Adapted from Pyomo.DAE: Racing Example Revisited and examples
Solver transcripts, mesh evidence, and reproducibility¶
The two tee=True outputs below are the full IPOPT transcripts for this example. They are useful when you need to inspect the transcribed NLP rather than only the ODE: with Radau collocation (ncp=3) and piecewise-constant controls, the coarse N=10 model has 153 variables and 144 equality constraints; the fine N=100 model has 1503 and 1404. The observed final times are 16.5783 s and 16.4711 s, respectively.
The conclusion is about the mesh, not a universal iteration-count rule: the coarse mesh traps the bang-bang switch inside one element, while the finer mesh reduces the transition element width. A moving-element formulation below aligns a full bang-bang switch with a boundary. The stored outputs show one successful run; rerunning the cells is the authoritative way to obtain solver output for your environment. The archived race-car-control and race-car-states result files let the figures be redrawn without rerunning IPOPT.
Read the physical model and time scaling, mesh comparison, moving-element equations, continuous reference, and initialization experiment. The full collocation derivation is in Dynamic Optimization with Collocation and Pyomo.DAE.
Learning Objectives¶
Introduces syntax for Pyomo.DAE
Shows time-scaling modeling trick
Practice using the units feature in Pyomo, and learn where in the workflow to check them (see Units in Pyomo.dae Models)
Install Packages and Load Modules¶
# Imports
import sys
if "google.colab" in sys.modules:
!wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
import helper
helper.easy_install()
else:
sys.path.insert(0, "../")
import helper
helper.set_plotting_style()import pyomo.environ as pyo
import matplotlib.pyplot as plt
import pyomo.dae as dae
from pyomo.environ import units
from pyomo.util.check_units import assert_units_consistentProblem Statement and Optimal Control Formulation¶
You are a race car driver with a simple goal. Drive distance in the minimal amount of time but come to a complete stop at the finish line.
Mathematically, you want to solve the following optimal control problem:
where is the acceleration/braking (your control variable) and is the drag coefficient (parameter).
For the numerical runs, m and m; is the applied acceleration in m/s, so has the same units. The code bounds m and , and requires s. These bounds are part of the stated model, not implicit consequences of time normalization.
Scale Time¶
Let where . Thus . The optimal control problem becomes:
This change of variables assumes . For a nonzero fixed start time use and multiply by the duration . Any explicit physical-time argument in the model must change too.
Declaring our Model with Pyomo.DAE¶
We can use Pyomo.dae to automatically formulate the collocation equations (i.e., add constraints that numerically integrate the ODE model).
def create_model1():
# Define the model
m = pyo.ConcreteModel()
# Define the model parameters
m.R = pyo.Param(initialize=0.001, units=1 / units.m) # Friction factor
m.L = pyo.Param(initialize=100.0, units=units.m) # Final position
# Define time
m.tau = dae.ContinuousSet(bounds=(0, 1)) # Dimensionless time set
m.tf = pyo.Var(initialize=1, bounds=(1e-3, None), units=units.s) # Final time
# Define remaining algebraic variables
m.x = pyo.Var(m.tau, bounds=(0, m.L + 50 * units.m), units=units.m) # Position
m.v = pyo.Var(m.tau, bounds=(0, None), units=units.m / units.s) # Velocity
m.u = pyo.Var(
m.tau, bounds=(-3.0, 1.0), initialize=0, units=units.m / units.s / units.s
) # Acceleration
# Define derivative variables.
#
# Units MUST be given explicitly: DerivativeVar defaults to dimensionless and
# does NOT inherit units from the state variable. Because tau is dimensionless,
# d/dtau has the SAME units as the state itself.
m.dx = dae.DerivativeVar(m.x, units=units.m)
m.dv = dae.DerivativeVar(m.v, units=units.m / units.s)
# Declare the objective (minimize final time)
m.obj = pyo.Objective(expr=m.tf)
# Define the constraints
# position
@m.Constraint(m.tau)
def ode1(m, i):
if i == 0:
return pyo.Constraint.Skip
return m.dx[i] == m.tf * m.v[i]
# velocity
@m.Constraint(m.tau)
def ode2(m, i):
if i == 0:
return pyo.Constraint.Skip
return m.dv[i] == m.tf * (m.u[i] - m.R * m.v[i] ** 2)
# Define the initial/boundary conditions
def _init(m):
yield m.x[0] == 0
yield m.x[1] == m.L
yield m.v[0] == 0
yield m.v[1] == 0
m.initcon = pyo.ConstraintList(rule=_init)
# Check dimensional consistency HERE, on the continuous model, BEFORE the
# discretization equations are added. See ./units_and_pyomo_dae.md and
# https://github.com/Pyomo/pyomo/issues/1790
assert_units_consistent(m)
return m
m = create_model1()Now let’s inspect the model.
m.pprint()2 Param Declarations
L : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=m
Key : Value
None : 100.0
R : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=1/m
Key : Value
None : 0.001
4 Var Declarations
tf : Size=1, Index=None, Units=s
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0.001 : 1 : None : False : False : Reals
u : Size=2, Index=tau, Units=m/s**2
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : -3.0 : 0 : 1.0 : False : False : Reals
1 : -3.0 : 0 : 1.0 : False : False : Reals
v : Size=2, Index=tau, Units=m/s
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : 0 : None : None : False : True : Reals
1 : 0 : None : None : False : True : Reals
x : Size=2, Index=tau, Units=m
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : 0 : None : 150.0 : False : True : Reals
1 : 0 : None : 150.0 : False : True : Reals
1 Objective Declarations
obj : Size=1, Index=None, Active=True
Key : Active : Sense : Expression
None : True : minimize : tf
3 Constraint Declarations
initcon : Size=4, Index={1, 2, 3, 4}, Active=True
Key : Lower : Body : Upper : Active
1 : 0.0 : x[0] : 0.0 : True
2 : L : x[1] : L : True
3 : 0.0 : v[0] : 0.0 : True
4 : 0.0 : v[1] : 0.0 : True
ode1 : Size=1, Index=tau, Active=True
Key : Lower : Body : Upper : Active
1 : 0.0 : dx[1] - tf*v[1] : 0.0 : True
ode2 : Size=1, Index=tau, Active=True
Key : Lower : Body : Upper : Active
1 : 0.0 : dv[1] - tf*(u[1] - R*v[1]**2) : 0.0 : True
1 ContinuousSet Declarations
tau : Size=1, Index=None, Ordered=Sorted
Key : Dimen : Domain : Size : Members
None : 1 : [0..1] : 2 : {0, 1}
2 DerivativeVar Declarations
dv : Size=2, Index=tau, Units=m/s
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals
dx : Size=2, Index=tau, Units=m
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals
13 Declarations: R L tau tf x v u dx dv obj ode1 ode2 initcon
Aside: checking units, and when to check them¶
Notice where assert_units_consistent(m) sits in create_model1 above: at the very
end of the model-building function, before any discretization is applied.
Two details in the code above are worth pausing on.
DerivativeVardoes not inherit units. It defaults to dimensionless. A model can look thoroughly annotated and still fail the very first check. Before this notebook declaredunits=onm.dxandm.dv,assert_units_consistentfailed immediately withdimensionless not compatible with meter-- even though everyVarandParamhad units.The units of are metres, not metres per second. We scaled time, so is dimensionless and the derivative carries the same units as the state.
That second point has a pleasant consequence, explored in the cell below.
See Units in Pyomo.dae Models for the full story, including what the check does
not look at (Var bounds and .fix() values are never checked).
# The time-scaling trick has a side benefit: because tau is DIMENSIONLESS, the
# finite-difference/collocation coefficients really are dimensionless, so this model
# passes the units check BEFORE and AFTER discretization.
#
# Do not generalize from this! It is a happy accident of scaling time to [0, 1].
# A model whose ContinuousSet is time-in-seconds fails after discretization --
# see notebooks/3/PyomoDAE_TCLab.ipynb for that case.
m_demo = create_model1() # assert_units_consistent already ran inside
print("Before discretization: units are consistent.")
pyo.TransformationFactory("dae.collocation").apply_to(
m_demo, nfe=15, scheme="LAGRANGE-RADAU", ncp=3
)
try:
assert_units_consistent(m_demo)
print("After discretization: units are STILL consistent (tau is dimensionless).")
except Exception as e:
print("After discretization: FAILED ->", e)Before discretization: units are consistent.
After discretization: units are STILL consistent (tau is dimensionless).
Discretize/Transcribe and Solve¶
Pyomo.dae writes the collocation equations for us. This experiment varies nfe, the number of finite elements, while holding three Radau stages per element and one control per element fixed. Changing stage count, control parameterization or mesh locations is a separate experiment. We solve the same model on
two meshes and compare.
Notice what tee=True prints before the iteration table -- Ipopt reports how
many variables and equality constraints the collocation equations produced.
That count, not the size of the ODE, is what the solver actually sees.
# COARSE mesh: N = 10 finite elements.
#
# N = 10 is chosen, not arbitrary. The optimal control here is bang-bang: full
# throttle (u = +1) until one switching time, then full braking (u = -3) to the
# line. At N = 10 that switching time falls strictly INSIDE a finite element,
# and because reduce_collocation_points holds u constant across an element, the
# only way the discretized problem can say "switch part-way through this
# element" is to give the WHOLE element an intermediate value. Watch for it
# below: one element sits at u = -0.52 m/s^2 for 1.7 s -- a braking level the
# true optimal control never uses. That is an artifact of the mesh, not physics.
NFE_COARSE = 10
m_coarse = create_model1()
# Declare the discretizer
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m_coarse, nfe=NFE_COARSE, scheme="LAGRANGE-RADAU", ncp=3)
# force piecewise constant controls (acceleration) over each finite element
m_coarse = discretizer.reduce_collocation_points(
m_coarse, var=m_coarse.u, ncp=1, contset=m_coarse.tau
)
# Solve. tee=True prints Ipopt's own log: first the problem-size block -- how
# many variables and equality constraints the collocation equations produced --
# and then the iteration table.
solver = pyo.SolverFactory("ipopt")
solve_status = solver.solve(m_coarse, tee=True)
assert pyo.check_optimal_termination(solve_status), (
f"Solve failed: status={solve_status.solver.status}, "
f"termination={solve_status.solver.termination_condition}"
)
print("final time = %6.2f seconds" % (pyo.value(m_coarse.tf)))Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 554
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 90
Total number of variables............................: 153
variables with only lower bounds: 32
variables with lower and upper bounds: 61
variables with only upper bounds: 0
Total number of equality constraints.................: 144
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 1.00e+02 2.79e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.5853340e+00 9.98e+01 1.13e+02 -1.0 1.52e+03 - 8.94e-05 1.70e-03f 1
2 3.7931592e+00 9.98e+01 1.12e+02 -1.0 2.32e+03 - 5.97e-04 4.89e-04h 2
3 4.0647486e+00 9.97e+01 1.12e+02 -1.0 2.18e+03 - 1.23e-03 6.67e-04h 1
4 4.1650879e+00 9.97e+01 1.11e+02 -1.0 1.23e+03 - 1.36e-03 3.14e-04h 1
5 4.6776122e+00 9.95e+01 1.09e+02 -1.0 1.12e+03 - 3.41e-03 1.65e-03h 1
6 5.2835145e+00 9.93e+01 1.07e+02 -1.0 8.83e+02 - 5.35e-03 2.19e-03h 1
7 6.0098618e+00 9.90e+01 1.04e+02 -1.0 7.24e+02 - 1.10e-04 2.95e-03h 2
8 6.3312839e+00 9.89e+01 1.04e+02 -1.0 6.59e+02 - 6.85e-03 1.46e-03h 1
9 7.3012932e+00 9.84e+01 1.00e+02 -1.0 5.84e+02 - 7.17e-03 4.72e-03h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 8.0540486e+00 9.80e+01 9.85e+01 -1.0 5.92e+02 - 2.69e-04 4.18e-03h 2
11 8.5947521e+00 9.77e+01 9.76e+01 -1.0 6.11e+02 - 5.37e-03 3.24e-03h 1
12 9.5877553e+00 9.70e+01 9.53e+01 -1.0 5.97e+02 - 2.95e-02 6.40e-03h 3
13 1.0448493e+01 9.64e+01 9.38e+01 -1.0 5.78e+02 - 6.44e-02 6.18e-03h 3
14 1.1589405e+01 9.56e+01 9.16e+01 -1.0 5.68e+02 - 5.49e-02 8.83e-03h 2
15 1.2853151e+01 9.46e+01 1.24e+02 -1.0 5.72e+02 - 3.05e-02 1.05e-02h 1
16 1.4009896e+01 9.36e+01 1.03e+02 -1.0 5.77e+02 - 1.58e-01 1.04e-02h 4
17 1.5265179e+01 9.25e+01 1.02e+02 -1.0 5.81e+02 - 1.42e-02 1.19e-02h 4
18 1.6627412e+01 9.12e+01 1.65e+02 -1.0 5.86e+02 - 2.42e-01 1.37e-02h 4
19 1.7931918e+01 9.00e+01 1.73e+02 -1.0 5.90e+02 - 6.13e-02 1.38e-02h 4
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.9226411e+01 8.87e+01 1.79e+02 -1.0 5.97e+02 - 5.13e-02 1.42e-02h 4
21 2.1876814e+01 8.60e+01 2.12e+02 -1.0 6.06e+02 - 1.27e-01 3.00e-02h 3
22 2.4847549e+01 8.30e+01 2.02e+02 -1.0 6.30e+02 - 3.63e-02 3.49e-02h 2
23 2.7281885e+01 8.05e+01 1.92e+02 -1.0 6.48e+02 - 2.56e-02 3.02e-02h 1
24 2.8056514e+01 7.93e+01 1.98e+02 -1.0 5.36e+02 - 1.03e-01 1.45e-02h 1
25 3.0856978e+01 7.49e+01 1.86e+02 -1.0 5.25e+02 - 8.92e-02 5.55e-02h 5
26 3.3534388e+01 7.03e+01 1.68e+02 -1.0 4.96e+02 - 1.83e-01 6.19e-02h 5
27 3.5839592e+01 6.59e+01 1.17e+02 -1.0 4.65e+02 - 4.60e-01 6.19e-02h 5
28 3.9870395e+01 5.78e+01 9.06e+01 -1.0 4.36e+02 - 3.62e-01 1.24e-01h 4
29 5.9406704e+01 1.43e+01 7.65e+01 -1.0 3.82e+02 - 9.98e-01 9.90e-01H 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 5.9207591e+01 1.87e-01 5.99e+02 -1.0 6.53e+00 - 6.44e-01 9.90e-01h 1
31 5.9028153e+01 5.92e-03 6.24e+04 -1.0 3.61e+00 - 6.39e-01 1.00e+00f 1
32 5.8294881e+01 7.32e-02 6.23e+06 -1.0 1.29e+01 - 3.70e-01 1.00e+00f 1
33 5.4933357e+01 1.68e+00 5.19e+06 -1.0 4.48e+01 - 1.68e-01 1.00e+00f 1
34 4.0577253e+01 1.74e+01 4.51e+06 -1.0 2.11e+02 - 1.31e-01 8.79e-01F 1
35 4.4448339e+01 6.40e+00 3.04e+06 -1.0 3.40e+01 - 3.25e-01 1.00e+00h 1
36 3.5629558e+01 6.10e+00 2.39e+06 -1.0 1.56e+02 - 2.16e-01 7.64e-01F 1
37 3.6657760e+01 1.14e+00 9.83e+05 -1.0 3.62e+01 - 5.88e-01 1.00e+00f 1
38 3.0907089e+01 3.72e+00 7.06e+05 -1.0 2.46e+02 - 2.81e-01 4.87e-01F 1
39 2.2125037e+01 3.11e+01 4.40e+05 -1.0 2.32e+02 - 3.77e-01 6.32e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 2.0255631e+01 2.21e+01 9.52e+04 -1.0 1.35e+02 - 7.84e-01 3.71e-01h 1
41 1.7838350e+01 1.31e+01 9.05e+03 -1.0 1.14e+02 - 9.05e-01 8.56e-01h 1
42 1.8616120e+01 9.01e-01 1.89e-02 -1.0 2.52e+01 - 1.00e+00 1.00e+00f 1
43 1.7357574e+01 1.63e+00 8.88e+05 -1.7 1.29e+01 - 8.89e-01 1.00e+00h 1
44 1.7079996e+01 1.36e-01 1.33e-03 -1.7 4.13e+00 - 1.00e+00 1.00e+00h 1
45 1.6685174e+01 2.01e-01 4.82e-03 -2.5 4.44e+00 - 1.00e+00 1.00e+00h 1
46 1.6653207e+01 2.77e-03 1.30e-05 -2.5 1.15e+00 - 1.00e+00 1.00e+00h 1
47 1.6583273e+01 6.65e-03 1.38e-04 -3.8 8.11e-01 - 1.00e+00 1.00e+00h 1
48 1.6582355e+01 1.88e-06 9.35e-09 -3.8 2.33e-02 - 1.00e+00 1.00e+00h 1
49 1.6578350e+01 2.24e-05 4.43e-07 -5.7 4.83e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 1.6578347e+01 1.82e-11 1.85e-11 -5.7 6.54e-05 - 1.00e+00 1.00e+00h 1
51 1.6578297e+01 3.46e-09 6.83e-11 -8.6 6.00e-04 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 51
(scaled) (unscaled)
Objective...............: 1.6578296960731734e+01 1.6578296960731734e+01
Dual infeasibility......: 6.8294368101657931e-11 6.8294368101657931e-11
Constraint violation....: 3.4584672903292812e-09 3.4584672903292812e-09
Complementarity.........: 2.5546384648403259e-09 2.5546384648403259e-09
Overall NLP error.......: 3.4584672903292812e-09 3.4584672903292812e-09
Number of objective function evaluations = 127
Number of objective gradient evaluations = 52
Number of equality constraint evaluations = 127
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 52
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 51
Total CPU secs in IPOPT (w/o function evaluations) = 0.012
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
final time = 16.58 seconds
Finer Time Discretization¶
Does the intermediate braking force go away if we consider a finer time discretization?
# FINE mesh: N = 100 finite elements.
#
# The same problem on a mesh ten times finer. Nothing about the model changes;
# only nfe does. The switching element is still there -- it always is, because
# the switch almost never lands exactly on an element boundary -- but it is now
# 0.16 s wide instead of 1.7 s, so the control profile looks like the bang-bang
# solution it is approximating, and the objective stops moving.
NFE_FINE = 100
m_fine = create_model1()
# Declare the discretizer
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m_fine, nfe=NFE_FINE, scheme="LAGRANGE-RADAU", ncp=3)
# force piecewise constant controls (acceleration) over each finite element
m_fine = discretizer.reduce_collocation_points(
m_fine, var=m_fine.u, ncp=1, contset=m_fine.tau
)
# Solve. tee=True prints Ipopt's own log: first the problem-size block -- how
# many variables and equality constraints the collocation equations produced --
# and then the iteration table.
solver = pyo.SolverFactory("ipopt")
solve_status = solver.solve(m_fine, tee=True)
assert pyo.check_optimal_termination(solve_status), (
f"Solve failed: status={solve_status.solver.status}, "
f"termination={solve_status.solver.termination_condition}"
)
print("final time = %6.2f seconds" % (pyo.value(m_fine.tf)))Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 5504
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 900
Total number of variables............................: 1503
variables with only lower bounds: 302
variables with lower and upper bounds: 601
variables with only upper bounds: 0
Total number of equality constraints.................: 1404
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 1.00e+02 1.70e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 8.3177905e+00 9.99e+01 5.12e+01 -1.0 1.05e+04 - 9.46e-05 7.00e-04f 1
2 9.1154780e+00 9.99e+01 5.09e+01 -1.0 4.36e+04 - 1.83e-05 1.83e-05s 7
3 9.7866770e+00 9.99e+01 5.07e+01 -1.0 1.81e+04 - 4.94e-05 4.94e-05s 7
4 1.0582478e+01 9.99e+01 5.47e+01 -1.0 1.37e+04 - 1.67e-04 0.00e+00S 7
5 1.2791004e+01 9.99e+01 5.25e+01 -1.0 1.10e+04 - 3.56e-04 3.02e-04h 1
6 1.3678029e+01 9.99e+01 5.22e+01 -1.0 1.03e+04 - 8.86e-04 1.50e-04h 2
7 1.5126950e+01 9.98e+01 6.38e+01 -1.0 9.52e+03 - 1.57e-03 2.65e-04h 2
8 1.6078088e+01 9.98e+01 1.29e+02 -1.0 8.75e+03 - 3.06e-03 1.92e-04h 3
9 1.7141305e+01 9.98e+01 2.16e+02 -1.0 7.49e+03 - 4.30e-03 2.29e-04h 4
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.8372431e+01 9.98e+01 3.23e+02 -1.0 6.53e+03 - 5.72e-03 2.82e-04h 4
11 1.8970184e+01 9.97e+01 3.88e+02 -1.0 5.72e+03 - 3.72e-03 1.47e-04h 6
12 1.9926487e+01 9.97e+01 5.77e+02 -1.0 5.41e+03 - 1.11e-02 2.43e-04h 6
13 2.0218795e+01 9.97e+01 6.34e+02 -1.0 4.88e+03 - 3.55e-03 7.75e-05h 9
14 2.0566503e+01 9.97e+01 9.54e+02 -1.0 4.77e+03 - 1.99e-02 9.34e-05h 9
15 2.0887021e+01 9.97e+01 1.01e+03 -1.0 4.38e+03 - 3.40e-03 8.74e-05h 10
16 1.9749372e+02 4.70e+02 3.54e+03 -1.0 4.31e+03 - 1.92e-02 4.88e-02w 1
17 1.8398997e+02 3.14e+02 1.52e+05 -1.0 1.13e+03 0.0 2.38e-04 2.97e-01w 1
18 3.1258025e+02 3.31e+02 1.07e+05 -1.0 6.88e+02 0.4 2.79e-04 3.23e-01w 1
19 2.1231956e+01 9.97e+01 1.30e+03 -1.0 1.83e+02 0.9 1.92e-02 9.54e-05h 9
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 2.1703484e+01 9.97e+01 1.35e+03 -1.0 4.11e+03 - 3.69e-03 1.32e-04h 10
21 2.2194551e+01 9.97e+01 2.09e+03 -1.0 4.04e+03 - 5.17e-02 1.41e-04h 10
22 2.2504323e+01 9.97e+01 2.13e+03 -1.0 3.78e+03 - 3.59e-03 9.06e-05h 11
23 2.2822532e+01 9.96e+01 2.30e+03 -1.0 3.74e+03 - 1.30e-02 9.42e-05h 11
24 2.3149526e+01 9.96e+01 2.36e+03 -1.0 3.68e+03 - 5.04e-03 9.80e-05h 11
25 2.3485674e+01 9.96e+01 4.57e+03 -1.0 3.64e+03 - 1.77e-01 1.02e-04h 11
26 2.3737411e+01 9.96e+01 4.61e+03 -1.0 3.42e+03 - 3.81e-03 7.77e-05h 11
27 2.4010006e+01 9.96e+01 4.71e+03 -1.0 3.39e+03 - 9.09e-03 8.49e-05h 11
28 2.4303959e+01 9.96e+01 4.80e+03 -1.0 3.36e+03 - 9.25e-03 9.25e-05h 11
29 3.5553336e+02 1.47e+03 9.23e+03 -1.0 3.33e+03 - 3.13e-02 1.05e-01w 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 3.1117768e+02 1.23e+03 1.42e+05 -1.0 7.80e+02 0.4 1.14e-03 1.60e-01w 1
31 4.3014533e+02 4.85e+02 4.69e+05 -1.0 3.87e+02 0.8 6.37e-05 9.18e-01w 1
32 2.4627426e+01 9.96e+01 5.11e+03 -1.0 8.97e+01 1.2 3.13e-02 1.03e-04h 10
33 2.4963404e+01 9.96e+01 5.15e+03 -1.0 3.28e+03 - 4.69e-03 1.08e-04h 11
34 2.5348201e+01 9.96e+01 8.02e+03 -1.0 3.25e+03 - 3.05e-01 1.25e-04h 11
35 2.5593752e+01 9.96e+01 8.05e+03 -1.0 3.12e+03 - 4.38e-03 8.17e-05h 11
36 2.5853840e+01 9.95e+01 8.12e+03 -1.0 3.09e+03 - 8.68e-03 8.73e-05h 11
37 2.6129670e+01 9.95e+01 8.28e+03 -1.0 3.07e+03 - 1.85e-02 9.34e-05h 11
38 2.6421370e+01 9.95e+01 8.37e+03 -1.0 3.04e+03 - 1.10e-02 9.98e-05h 11
39 2.6733933e+01 9.95e+01 8.43e+03 -1.0 3.01e+03 - 8.11e-03 1.08e-04h 11
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 2.7073242e+01 9.95e+01 8.77e+03 -1.0 2.98e+03 - 4.39e-02 1.19e-04h 11
41 2.7424679e+01 9.95e+01 8.81e+03 -1.0 2.94e+03 - 5.51e-03 1.24e-04h 11
42 4.2508223e+02 2.10e+03 1.05e+04 -1.0 2.90e+03 - 1.81e-01 1.42e-01w 1
43 2.3386965e+02 5.61e+02 1.38e+06 -1.0 6.64e+02 0.7 9.27e-04 7.41e-01w 1
44 3.7855705e+02 6.96e+02 7.48e+05 -1.0 4.64e+02 1.2 2.54e-02 8.43e-01w 1
45 2.7813016e+01 9.95e+01 1.02e+04 -1.0 1.26e+02 2.5 1.81e-01 1.39e-04h 10
46 2.8158499e+01 9.95e+01 1.02e+04 -1.0 2.84e+03 - 5.31e-03 1.26e-04h 11
47 2.8533873e+01 9.95e+01 1.03e+04 -1.0 2.81e+03 - 1.25e-02 1.38e-04h 11
48 2.8943932e+01 9.94e+01 1.04e+04 -1.0 2.78e+03 - 1.14e-02 1.53e-04h 11
49 2.9397995e+01 9.94e+01 1.06e+04 -1.0 2.74e+03 - 2.40e-02 1.71e-04h 11
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 2.9649224e+01 9.94e+01 1.07e+04 -1.0 2.70e+03 - 7.51e-03 9.63e-05h 12
51 2.9917372e+01 9.94e+01 1.46e+04 -1.0 2.68e+03 - 4.92e-01 1.04e-04h 12
52 3.0326448e+01 9.94e+01 1.47e+04 -1.0 2.62e+03 - 6.00e-03 1.60e-04h 11
53 3.0768578e+01 9.94e+01 1.50e+04 -1.0 2.59e+03 - 1.65e-02 1.75e-04h 11
54 3.1010574e+01 9.94e+01 1.54e+04 -1.0 2.55e+03 - 2.49e-02 9.74e-05h 12
55 5.5461600e+02 3.67e+03 2.25e+04 -1.0 2.53e+03 - 1.31e-02 2.12e-01w 1
56 2.6548248e+02 1.13e+03 1.16e+07 -1.0 3.80e+02 2.0 2.93e-04 7.61e-01w 1
57 1.9780854e+02 1.66e+02 2.19e+06 -1.0 1.53e+02 2.5 5.44e-02 9.90e-01h 1
58 1.6962504e+02 8.72e+00 7.25e+05 -1.0 2.85e+01 2.9 6.94e-01 9.91e-01h 1
59 1.6892261e+02 3.88e-02 3.02e+04 -1.0 1.84e+00 2.4 9.44e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 1.6891830e+02 1.32e-06 1.35e+00 -1.0 5.52e-03 1.9 1.00e+00 1.00e+00h 1
61 1.6891567e+02 1.55e-07 1.22e-01 -3.8 4.36e-03 1.4 1.00e+00 1.00e+00h 1
62 1.6890777e+02 1.39e-06 1.22e-01 -3.8 1.31e-02 1.0 1.00e+00 1.00e+00f 1
63 1.6888407e+02 1.25e-05 1.22e-01 -3.8 3.92e-02 0.5 1.00e+00 1.00e+00f 1
64 1.6881305e+02 1.13e-04 1.22e-01 -3.8 1.18e-01 0.0 1.00e+00 1.00e+00f 1
65 1.6860056e+02 1.01e-03 1.22e-01 -3.8 3.52e-01 -0.5 1.00e+00 1.00e+00f 1
66 1.6796823e+02 8.96e-03 1.21e-01 -3.8 1.05e+00 -0.9 1.00e+00 1.00e+00f 1
67 1.6611621e+02 7.77e-02 1.19e-01 -3.8 3.09e+00 -1.4 1.00e+00 1.00e+00f 1
68 1.6093290e+02 6.28e-01 1.13e-01 -3.8 8.77e+00 -1.9 1.00e+00 1.00e+00f 1
69 1.4805568e+02 4.23e+00 1.34e-01 -3.8 2.27e+01 -2.4 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 1.3617126e+02 6.75e+00 1.75e-01 -3.8 4.88e+01 -2.8 1.00e+00 4.79e-01f 1
71 1.3616854e+02 6.75e+00 6.87e+00 -3.8 6.01e+01 -3.3 1.00e+00 2.54e-04h 1
72 1.1510892e+02 1.80e+01 1.68e+00 -3.8 1.88e+02 -3.8 1.00e+00 1.00e+00f 1
73 8.8758484e+01 4.81e+01 1.86e+00 -3.8 5.56e+02 -4.3 6.15e-01 3.20e-01f 1
74 8.7921597e+01 4.80e+01 4.03e+00 -3.8 3.17e+03 -4.8 7.81e-02 1.87e-03h 1
75 8.4984976e+01 4.68e+01 4.51e+00 -3.8 3.60e+02 -4.3 2.86e-01 3.99e-02h 1
76 8.3941971e+01 4.68e+01 4.46e+00 -3.8 3.56e+04 -4.8 8.14e-07 3.51e-04f 1
77 8.0907192e+01 4.28e+01 4.08e+00 -3.8 4.18e+02 -4.4 1.00e+00 8.51e-02f 1
78 8.0822358e+01 4.24e+01 3.06e+00 -3.8 1.22e+02 -4.0 1.00e+00 8.84e-03h 1
79 6.2491310e+01 4.40e+01 2.22e+00 -3.8 5.02e+02 -4.4 5.00e-01 4.09e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 6.2432089e+01 4.40e+01 2.54e+00 -3.8 6.53e+03 -4.9 1.47e-02 1.08e-04h 1
81 5.2917043e+01 3.26e+01 1.94e+00 -3.8 4.84e+02 -4.5 8.75e-02 3.29e-01h 1
82 5.2914079e+01 3.26e+01 1.92e+00 -3.8 8.98e+01 -4.1 5.64e-01 4.50e-04h 1
83 3.1746112e+01 8.19e+01 7.37e-01 -3.8 3.77e+02 -4.5 1.94e-03 8.79e-01h 1
84 3.1744026e+01 8.19e+01 6.20e-01 -3.8 1.13e+02 -4.1 3.31e-01 3.20e-04h 1
85 3.0109091e+01 7.74e+01 6.41e-01 -3.8 3.21e+02 -4.6 2.44e-01 5.55e-02h 1
86 2.9642779e+01 7.23e+01 6.00e-01 -3.8 9.91e+01 -4.2 1.62e-01 6.60e-02h 1
87 2.7500336e+01 6.47e+01 5.40e-01 -3.8 3.12e+02 -4.6 1.05e-01 1.04e-01h 1
88 2.6879769e+01 5.56e+01 4.66e-01 -3.8 8.14e+01 -4.2 1.00e+00 1.41e-01h 1
89 2.4712134e+01 5.12e+01 4.31e-01 -3.8 3.57e+02 -4.7 9.20e-02 8.02e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 2.4177948e+01 4.61e+01 3.90e-01 -3.8 7.61e+01 -4.3 1.00e+00 9.85e-02h 1
91 2.1877168e+01 4.20e+01 3.57e-01 -3.8 3.24e+02 -4.7 1.14e-01 8.91e-02h 1
92 2.1818916e+01 4.13e+01 3.51e-01 -3.8 6.89e+01 -4.3 2.76e-01 1.60e-02h 1
93 2.1663311e+01 4.11e+01 3.49e-01 -3.8 4.18e+02 -4.8 1.59e-01 6.40e-03h 1
94 1.9268982e+01 4.03e+01 3.38e-01 -3.8 1.73e+03 -5.3 5.78e-03 1.89e-02h 1
95 1.9021233e+01 3.98e+01 3.34e-01 -3.8 5.28e+02 -4.8 4.13e-02 1.20e-02h 1
96 1.8994105e+01 3.92e+01 3.50e-01 -3.8 7.08e+01 -4.4 2.79e-01 1.48e-02h 1
97 1.8286769e+01 3.70e+01 3.34e-01 -3.8 1.96e+02 -4.9 3.47e-02 5.53e-02h 1
98 1.8191455e+01 3.49e+01 4.03e+00 -3.8 6.30e+01 -4.5 1.43e-03 5.86e-02h 1
99 1.7915228e+01 3.40e+01 3.88e+00 -3.8 2.64e+02 -4.9 1.02e-01 2.40e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 1.7758543e+01 2.94e+01 3.44e+00 -3.8 5.69e+01 -4.5 2.70e-01 1.36e-01h 1
101 1.7320249e+01 2.69e+01 3.12e+00 -3.8 1.19e+02 -5.0 2.78e-01 8.58e-02h 1
102 1.7243554e+01 2.39e+01 2.78e+00 -3.8 6.19e+01 -4.6 1.31e-01 1.11e-01h 1
103 1.6915068e+01 2.21e+01 2.57e+00 -3.8 1.06e+02 -5.0 6.60e-02 7.48e-02h 1
104 1.6874609e+01 2.10e+01 2.44e+00 -3.8 5.55e+01 -4.6 2.94e-02 4.95e-02h 1
105 1.6576222e+01 1.95e+01 2.27e+00 -3.8 1.05e+02 -5.1 2.04e-01 7.16e-02h 1
106 1.6397977e+01 1.52e+01 1.77e+00 -3.8 6.73e+01 -4.7 1.13e-01 2.22e-01h 1
107 1.6193458e+01 1.41e+01 1.64e+00 -3.8 7.54e+01 -5.1 2.73e-01 7.10e-02h 1
108 1.6150903e+01 1.41e+01 1.64e+00 -3.8 4.47e+04 - 4.60e-04 2.08e-04h 1
109 1.6156439e+01 1.40e+01 1.63e+00 -3.8 7.72e+03 - 1.34e-02 7.47e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 1.6053857e+01 1.26e+01 1.47e+00 -3.8 2.58e+02 - 7.69e-02 9.86e-02h 1
111 1.6069574e+01 1.19e+01 1.39e+00 -3.8 5.44e+02 - 3.67e-02 5.43e-02h 1
112 1.6056315e+01 1.09e+01 1.27e+00 -3.8 1.60e+02 - 4.56e-01 8.51e-02h 1
113 1.6255365e+01 5.49e+00 6.38e-01 -3.8 4.50e+01 - 1.00e+00 4.97e-01h 1
114 1.6515746e+01 3.95e-02 1.06e-03 -3.8 5.91e+00 - 1.00e+00 1.00e+00h 1
115 1.6515472e+01 1.39e-06 2.11e-07 -3.8 1.23e-01 - 1.00e+00 1.00e+00h 1
116 1.6472035e+01 1.31e-02 1.28e+03 -5.7 4.96e+00 - 9.14e-01 1.00e+00f 1
117 1.6471688e+01 6.86e-06 2.64e-09 -5.7 3.28e-01 - 1.00e+00 1.00e+00h 1
118 1.6471138e+01 1.80e-06 7.87e-09 -8.6 5.26e-02 - 1.00e+00 1.00e+00h 1
119 1.6471138e+01 1.07e-11 2.55e-14 -8.6 1.18e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 119
(scaled) (unscaled)
Objective...............: 1.6471137966472494e+01 1.6471137966472494e+01
Dual infeasibility......: 2.5545174048500027e-14 2.5545174048500027e-14
Constraint violation....: 1.8782218919210525e-12 1.0662682389876868e-11
Complementarity.........: 2.5059126720373543e-09 2.5059126720373543e-09
Overall NLP error.......: 2.5059126720373543e-09 2.5059126720373543e-09
Number of objective function evaluations = 538
Number of objective gradient evaluations = 120
Number of equality constraint evaluations = 538
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 120
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 119
Total CPU secs in IPOPT (w/o function evaluations) = 0.316
Total CPU secs in NLP function evaluations = 0.025
EXIT: Optimal Solution Found.
final time = 16.47 seconds
Comparing the Two Meshes¶
Both solves are done, so pull the trajectories out of the two models and
archive them. Separating “solve” from “analyze” is a habit worth keeping: the
plotting cells below then depend only on results, and can be re-run without a
solver.
def extract_race_car_results(m, nfe):
"""Pull the solved trajectories out of a Pyomo model into plain lists.
Arguments:
m: solved, discretized race-car model
nfe: number of finite elements used, for the label
Returns:
dict of JSON-safe values -- everything the figures need and nothing else
The time axis is rebuilt here: the model is written in dimensionless time
tau in [0, 1], so physical time is t = tau * t_f, and t_f is itself a
decision variable. Two meshes therefore have DIFFERENT final times, which
is the number the mesh-refinement argument turns on.
"""
tf = pyo.value(m.tf)
tau = sorted(m.tau)
return {
"label": f"N = {nfe}",
"nfe": nfe,
"tf": tf,
"t": [i * tf for i in tau], # time [s]
"x": [pyo.value(m.x[i]) for i in tau], # position [m]
"v": [pyo.value(m.v[i]) for i in tau], # velocity [m/s]
"u": [pyo.value(m.u[i]) for i in tau], # acceleration [m/s^2]
# Element boundaries in physical time, so a figure can show WHERE the
# control is allowed to change value.
"elements": [b * tf for b in m.tau.get_finite_elements()],
"element_u": [pyo.value(m.u[b]) for b in m.tau.get_finite_elements()[1:]],
}
results = {
"meshes": [
extract_race_car_results(m_coarse, NFE_COARSE),
extract_race_car_results(m_fine, NFE_FINE),
],
"u_lower": -3.0, # [m/s^2]
"u_upper": 1.0, # [m/s^2]
"L": 100.0, # [m]
}
for mesh in results["meshes"]:
print(f"{mesh['label']:>8}: t_f = {mesh['tf']:.4f} s")
DESCRIPTION = (
"Minimum-time race car (m = 1, L = 100 m, drag R = 0.001 1/m, "
"-3 <= u <= 1), solved with Pyomo.dae LAGRANGE-RADAU collocation, ncp = 3, "
"with the control reduced to piecewise constant on each element. Two "
"meshes: N = 10, coarse enough that the bang-bang switch is trapped inside "
"one element and that element takes an intermediate braking value, and "
"N = 100, where the switch is resolved."
)
for figure_name in ("race-car-control", "race-car-states"):
helper.save_results(
figure_name,
results,
notebook="notebooks/3-dev/PyomoDAE_car.ipynb",
source_tag="handout:race-car-model",
description=DESCRIPTION,
solver="Ipopt 3.13.2 (IDAES build, linear solver ma27)",
source_cells=[("notebooks/3-dev/PyomoDAE_car.ipynb", tag) for tag in
("handout-output:race-car-coarse", "handout-output:race-car-fine",
"numeric:race-car-extract")],
) N = 10: t_f = 16.5783 s
N = 100: t_f = 16.4711 s
[helper] wrote figures/results/race-car-control.json
[helper] wrote figures/results/race-car-states.json
The control profile¶
This is where the two meshes disagree.
# This cell is tagged `figure:race-car-control', which makes it the single
# source of the control-profile figure the Lecture 7 handout prints.
# figures/render_from_notebook.py re-runs exactly this cell against the archived
# JSON when the house style changes, so there is no second copy of the plotting
# code. That driver binds `np', `plt' and `helper' for the cell, so it imports
# nothing and uses no name defined in an earlier cell except `results'.
def plot_race_car_control(results):
"""The control profile on both meshes, and what the coarse mesh gets wrong.
THE FIGURE'S POINT is the flat step at u ~ -0.5 in the coarse profile. The
exact solution is bang-bang -- it takes only the two bound values -- so any
element sitting between them is the discretization confessing that it put
the switching time inside an element and had to average across it. On the
fine mesh the same element is a tenth as wide and the profile reads as the
step function it approximates.
STEPS AT ELEMENT EDGES. The control is piecewise constant on each finite
element. stairs(values, edges) places every jump at the element boundary;
stepping through collocation times would shift a jump into the next element.
⚠ tau = 0 IS DROPPED from the control series. u[0] is not a control the
solver uses -- the reduced control is indexed by element, and the first
element's value lives at its first collocation point -- so plotting u[0]
puts a spurious step from its initialize=0 value at the left edge.
ELEMENT RULES, coarse mesh only. The light vertical lines are the finite
element boundaries: the control is constant between two of them by
construction, so the rules explain the staircase rather than decorating it.
Drawing them for N = 100 would be a grey block.
GREYSCALE. Two series, each taking colour AND linestyle together from the
house prop_cycle (figures/dowling.mplstyle pairs the two element-wise),
plus different line widths, plus a direct label beside each curve. Any one
of those is enough on a photocopy.
"""
meshes = results["meshes"]
coarse = meshes[0]
u_lo, u_hi = results["u_lower"], results["u_upper"]
t_max = max(m["tf"] for m in meshes)
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
fig, ax = plt.subplots(figsize=(9.0, 4.2))
# Element boundaries FIRST, so they sit under the data.
for b in coarse["elements"]:
ax.axvline(b, color="0.85", linewidth=0.8, zorder=0)
for u_bound in (u_hi, u_lo):
ax.axhline(
u_bound, color="0.55", linewidth=1.0, linestyle=(0, (1, 3)), zorder=1
)
for k, mesh in enumerate(meshes):
# Radau stores an element's constant control at its right endpoint.
# The jump belongs at the shared element boundary, not at the next
# element's first interior collocation point.
ax.stairs(mesh["element_u"], mesh["elements"], baseline=None,
linewidth=3.2 if k == 0 else 1.8, zorder=3-k,
linestyle="-" if k == 0 else "--", color=colors[k])
# The defect, called out where it happens. Exactly one coarse element sits
# strictly between the bounds; find it rather than hard-coding a time.
tol = 1e-3
trapped = [
(t, u)
for t, u in zip(coarse["t"][1:], coarse["u"][1:])
if u_lo + tol < u < u_hi - tol
]
if trapped:
t_bad = trapped[len(trapped) // 2][0]
u_bad = trapped[0][1]
ax.annotate(
"the switch is trapped inside one element:\n"
f"the solver can only average across it, {u_bad:+.2f} m/s$^2$\n"
f"held for the element's full {coarse['tf'] / coarse['nfe']:.1f} s",
xy=(t_bad, u_bad),
xytext=(0.13, 0.30),
textcoords="axes fraction",
ha="left",
va="center",
fontsize=11,
arrowprops=dict(
arrowstyle="->",
linewidth=1.1,
color="0.25",
connectionstyle="arc3,rad=-0.25",
),
)
# DIRECT LABELS, not a legend: the two profiles coincide for most of the
# horizon, so a legend box would sit somewhere the eye has to travel to.
# They go ABOVE the engine-limit line -- headroom made by the ylim below,
# because every other corner of these axes has a curve or a rule in it.
for k, mesh in enumerate(meshes):
ax.annotate(
f"{mesh['label']}, $t_f$ = {mesh['tf']:.2f} s",
xy=(0.985, 0.95 - 0.09 * k),
xycoords="axes fraction",
ha="right",
va="center",
fontsize=12.5,
color=colors[k],
)
ax.annotate(
"engine limit, $u = +1$",
xy=(0.4, u_hi - 0.12),
ha="left",
va="top",
fontsize=10.5,
color="0.35",
)
ax.annotate(
"brake limit, $u = -3$",
xy=(0.4, u_lo + 0.14),
ha="left",
va="bottom",
fontsize=10.5,
color="0.35",
)
ax.annotate(
f"vertical rules: the {coarse['nfe']} finite element boundaries",
xy=(0.4, u_lo - 0.62),
ha="left",
va="bottom",
fontsize=10.5,
color="0.45",
)
ax.set_xlabel("time (s)")
ax.set_ylabel(r"acceleration $u$ (m/s$^2$)")
ax.set_xlim(0, t_max * 1.02)
ax.set_ylim(u_lo - 0.8, u_hi + 1.5)
fig.tight_layout()
return fig
fig = plot_race_car_control(results)
# Write media/figures/race-car-control.{png,pdf} -- what the lecture handout
# \includegraphics. A no-op on Colab, where there is no repo to write to.
helper.save_figure(fig, "race-car-control")[helper] wrote media/figures/race-car-control.png and .pdf
The state trajectories¶
And this is where they agree -- which is the trap.
# Tagged `figure:race-car-states': the single source of the state-trajectory
# figure in the Lecture 7 handout. Same contract as the cell above -- it uses
# only `results', `plt' and `helper', and defines its own figsize.
def plot_race_car_states(results):
"""Position and velocity on both meshes.
WHY THIS FIGURE SITS BESIDE THE CONTROL FIGURE. That one shows where the
coarse discretization goes wrong; this one shows how little of it reaches
the states. Position is visually identical on the two meshes, and the coarse
solution still lands exactly on x = L with v = 0 -- because those are
constraints, and a constraint holds on any mesh. What the coarse mesh gets
wrong is the OBJECTIVE. A student reading only these two panels would
conclude N = 10 was fine.
THE VELOCITY PANEL IS WHERE THE DEFECT IS VISIBLE. The trapped element's
weak braking rounds the peak of the coarse velocity curve where the fine one
turns sharply, and the coarse car then brakes at the limit for longer to
still stop on the line.
GREYSCALE. Colour and linestyle move together from the house prop_cycle,
the line widths differ, and the labels are direct rather than a legend.
"""
meshes = results["meshes"]
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
t_max = max(m["tf"] for m in meshes)
fig, (ax_x, ax_v) = plt.subplots(1, 2, figsize=(11.0, 4.0))
ax_x.axhline(
results["L"], color="0.55", linewidth=1.0, linestyle=(0, (1, 3)), zorder=1
)
for k, mesh in enumerate(meshes):
style = dict(linewidth=3.2 if k == 0 else 1.8, zorder=3 - k)
ax_x.plot(mesh["t"], mesh["x"], **style)
ax_v.plot(mesh["t"], mesh["v"], **style)
ax_x.annotate(
f"finish line, $L$ = {results['L']:.0f} m",
xy=(0.4, results["L"] - 3),
ha="left",
va="top",
fontsize=10.5,
color="0.35",
)
# Both labels go in the empty upper-left wedge of the velocity panel, above
# the acceleration ramp and left of the peak.
for k, mesh in enumerate(meshes):
ax_v.annotate(
f"{mesh['label']}, $t_f$ = {mesh['tf']:.2f} s",
xy=(0.03, 0.93 - 0.09 * k),
xycoords="axes fraction",
ha="left",
va="center",
fontsize=12.5,
color=colors[k],
)
ax_v.annotate(
"stopped on the line:\n$v(t_f) = 0$ is a constraint,\nso it holds on both meshes",
xy=(meshes[0]["tf"], 0.0),
xytext=(0.42, 0.30),
textcoords="axes fraction",
ha="left",
va="center",
fontsize=10.5,
color="0.35",
arrowprops=dict(
arrowstyle="->",
linewidth=1.0,
color="0.45",
connectionstyle="arc3,rad=-0.2",
),
)
ax_x.set_xlabel("time (s)")
ax_x.set_ylabel("position $x$ (m)")
ax_v.set_xlabel("time (s)")
ax_v.set_ylabel("velocity $v$ (m/s)")
for ax in (ax_x, ax_v):
ax.set_xlim(0, t_max * 1.02)
ax_v.set_ylim(0, max(max(m["v"]) for m in meshes) * 1.30)
fig.tight_layout()
return fig
fig = plot_race_car_states(results)
helper.save_figure(fig, "race-car-states")[helper] wrote media/figures/race-car-states.png and .pdf
Interpreting the mesh comparison¶
The transcribed endpoint constraints hold to the solver tolerance on both meshes. They do not certify the continuous trajectory or control-switch location. Compare final time, control values, the mesh, and a separately integrated trajectory. The moving-element experiment below supplies an independent bang-bang reference.
Initialization on a stated model¶
The bounded-duration model uses final time at least 0.001 seconds. The reproducible experiment below uses 15 elements, three Radau points, and controls free at every collocation point, so its results are not directly comparable with the piecewise-constant-control mesh runs. It records the initial duration, duration bound, final time and solver iterations. Iteration counts depend on the software and linear solver; a poor guess is not a guaranteed failure.
For harder problems, bound physical durations, initialize from simulation when possible, and solve coarse before refining.
The deliberately out-of-bounds zero-duration starting value tests how Ipopt projects an invalid guess into the stated positive-duration domain. It does not relax that bound. All rows use the same formulation.
An initial value that makes a logarithm, reciprocal, or other model expression undefined can prevent even the first residual evaluation. Stiffness and poor variable or residual scaling can also make the first local model unhelpful. The present experiment alone does not identify why recovery succeeds; calling the model “nearly linear” is not an explanation of every nonlinear solve. Fix controls at a nominal profile and use a successful simulation to initialize the states when possible. Simulator(m).initialize_model() supplies that pattern; the temperature-control lab example gives the full workflow. Check that the model and simulator support the intended DAE before relying on it.
Optimize ten finite-element lengths¶
The original model, coefficients, state bounds and control bounds are retained. This explicit Radau transcription makes physical element lengths decision variables, with nominal length 1.65 s and a 45 percent movement band: 0.9075 to 2.3925 s. Their sum is the final time. Position and velocity remain continuous; one control value belongs to each element.
For stage in element , use obtained by integrating the three-point Radau derivative basis:
The last Radau stage is the endpoint. Warm-start from the uniform ten-element solution. The control switch can then move onto a shared element boundary; this is a different formulation from changing only the nfe option in Pyomo.DAE.
Biegler, Section 10.3.1, equations (10.20)-(10.22), pp. 296-298.
More generally, Biegler’s movement band is , with positive nominal lengths and . These bounds prevent collapse and limit mesh movement. This experiment uses s and . For normalized fractions, use and instead of applying the duration twice. Fix the lengths for the first solve before releasing them.
Setting removes within-element control oscillations and reduces the control dimension. It is a modeling restriction, not a universal cure for singular controls: it can also prevent a switch inside an element. Moving boundaries and increasing the element count address that restriction in different ways.
import numpy as np
from numpy.polynomial import Polynomial
def radau_three_coefficients():
c = np.array([(4-np.sqrt(6))/10, (4+np.sqrt(6))/10, 1.0])
A = np.zeros((3, 3))
for k in range(3):
basis = Polynomial([1.0])
for j in range(3):
if j != k:
basis *= Polynomial([-c[j], 1.0]) / (c[k]-c[j])
integral = basis.integ()
A[:, k] = integral(c)-integral(0.0)
assert np.allclose(A.sum(axis=1), c)
assert np.isclose(A[-1].sum(), 1.0)
return c, A
def create_variable_element_car(nfe=10, nominal=1.65, gamma=0.45):
c, A = radau_three_coefficients()
m = pyo.ConcreteModel()
m.I = pyo.RangeSet(1, nfe)
m.J = pyo.RangeSet(0, 2)
m.B = pyo.RangeSet(0, nfe)
m.R = pyo.Param(initialize=0.001, units=1/units.m)
m.L = pyo.Param(initialize=100.0, units=units.m)
m.h = pyo.Var(m.I, bounds=((1-gamma)*nominal, (1+gamma)*nominal),
initialize=nominal, units=units.s)
m.tf = pyo.Var(bounds=(1e-3, None), initialize=nfe*nominal, units=units.s)
m.x = pyo.Var(m.I, m.J, bounds=(0, m.L+50*units.m), units=units.m)
m.v = pyo.Var(m.I, m.J, bounds=(0, None), units=units.m/units.s)
m.u = pyo.Var(m.I, bounds=(-3, 1), units=units.m/units.s**2)
m.xb = pyo.Var(m.B, bounds=(0, m.L+50*units.m), units=units.m)
m.vb = pyo.Var(m.B, bounds=(0, None), units=units.m/units.s)
m.xb[0].fix(0); m.vb[0].fix(0)
m.xb[nfe].fix(100); m.vb[nfe].fix(0)
@m.Constraint()
def horizon(m):
return sum(m.h[i] for i in m.I) == m.tf
@m.Constraint(m.I, m.J)
def position(m, i, j):
return m.x[i,j] == m.xb[i-1] + m.h[i]*sum(A[j,k]*m.v[i,k] for k in m.J)
@m.Constraint(m.I, m.J)
def velocity(m, i, j):
return m.v[i,j] == m.vb[i-1] + m.h[i]*sum(
A[j,k]*(m.u[i]-m.R*m.v[i,k]**2) for k in m.J)
@m.Constraint(m.I)
def position_end(m, i):
return m.xb[i] == m.x[i,2]
@m.Constraint(m.I)
def velocity_end(m, i):
return m.vb[i] == m.v[i,2]
m.obj = pyo.Objective(expr=m.tf)
assert_units_consistent(m)
# A feasible physical shape from the already solved uniform mesh.
old = extract_race_car_results(m_coarse, NFE_COARSE)
initial_edges = np.arange(nfe+1)*nominal
for i in m.B:
time = initial_edges[i]/initial_edges[-1]*old["tf"]
m.xb[i].set_value(np.interp(time, old["t"], old["x"]))
m.vb[i].set_value(np.interp(time, old["t"], old["v"]))
# Reset fixed endpoints exactly after interpolation.
m.xb[0].fix(0); m.vb[0].fix(0)
m.xb[nfe].fix(100); m.vb[nfe].fix(0)
for i in m.I:
m.u[i].set_value(old["element_u"][min(i-1,len(old["element_u"])-1)])
for j in m.J:
time = (initial_edges[i-1]+nominal*c[j])/initial_edges[-1]*old["tf"]
m.x[i,j].set_value(np.interp(time, old["t"], old["x"]))
m.v[i,j].set_value(np.interp(time, old["t"], old["v"]))
return m
m_variable = create_variable_element_car()
variable_status = pyo.SolverFactory("ipopt").solve(m_variable, tee=True,
options={"tol":1e-9,"max_iter":1000})
assert pyo.check_optimal_termination(variable_status)
print("Moving elements:", [pyo.value(m_variable.h[i]) for i in m_variable.I])
print("Controls:", [pyo.value(m_variable.u[i]) for i in m_variable.I])
print("Final time:", pyo.value(m_variable.tf))
Ipopt 3.13.2: tol=1e-09
max_iter=1000
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 403
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 70
Total number of variables............................: 99
variables with only lower bounds: 40
variables with lower and upper bounds: 59
variables with only upper bounds: 0
Total number of equality constraints.................: 81
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.6500000e+01 8.28e-02 1.86e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.7196867e+01 4.41e-02 5.51e-01 -1.0 7.04e-01 - 9.34e-01 9.90e-01f 1
2 1.7175451e+01 1.55e-01 1.70e+01 -1.0 1.92e+00 - 9.73e-01 9.90e-01f 1
3 1.7153115e+01 2.93e-01 3.32e+03 -1.0 1.63e+01 - 4.84e-01 5.17e-01f 1
4 1.7135748e+01 6.97e-02 6.65e-02 -1.0 4.27e+00 - 1.00e+00 1.00e+00h 1
5 1.6721031e+01 4.95e-02 3.50e+03 -2.5 7.21e-01 - 8.81e-01 1.00e+00h 1
6 1.6691228e+01 2.59e-01 3.22e+03 -2.5 8.65e+01 -4.0 7.99e-02 6.45e-02h 1
7 1.6517894e+01 1.27e-02 5.29e-03 -2.5 9.96e-01 - 1.00e+00 9.62e-01h 1
8 1.6472417e+01 6.88e-04 1.36e-04 -3.8 5.03e-01 - 1.00e+00 1.00e+00h 1
9 1.6470375e+01 7.21e-06 1.80e-06 -5.7 1.46e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.6470356e+01 5.71e-05 8.81e-01 -8.6 3.24e-02 - 9.94e-01 9.99e-01h 1
11 1.6470356e+01 1.57e-01 7.22e-01 -8.6 3.63e+01 - 1.53e-01 4.69e-02h 1
12r 1.6470356e+01 1.57e-01 9.99e+02 -0.8 0.00e+00 -2.7 0.00e+00 9.46e-10R 7
13r 1.6626010e+01 1.86e-01 9.99e+02 -0.8 7.50e+06 - 1.39e-12 3.40e-07f 1
14r 1.6615282e+01 1.32e-01 9.97e+02 -0.8 2.55e+02 - 2.39e-03 9.97e-04f 1
15r 1.6645671e+01 1.13e-01 9.91e+02 -0.8 1.08e+02 - 5.62e-03 3.52e-03f 1
16r 1.6655686e+01 1.02e-01 7.25e+02 -0.8 1.08e+01 - 2.69e-01 1.03e-02f 1
17r 1.6798281e+01 8.67e-02 5.94e+02 -0.8 2.34e+01 - 1.81e-01 4.04e-02f 1
18r 1.7005737e+01 4.85e-02 7.53e+01 -0.8 2.73e+00 - 1.00e+00 4.25e-01f 1
19r 1.7420285e+01 3.24e-02 4.48e-02 -0.8 4.78e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 1.7127839e+01 2.45e-02 5.43e-01 -3.4 5.33e-01 - 8.68e-01 1.00e+00f 1
21r 1.7114344e+01 1.27e-04 2.90e-03 -3.4 8.11e-02 - 1.00e+00 1.00e+00h 1
22r 1.6826497e+01 1.40e-02 8.65e-01 -5.1 3.68e-01 - 6.75e-01 1.00e+00f 1
23r 1.6803381e+01 1.18e-04 6.69e-05 -5.1 3.72e-02 - 1.00e+00 1.00e+00h 1
24r 1.6662449e+01 2.30e-03 1.18e-05 -5.1 1.68e-01 - 1.00e+00 1.00e+00h 1
25r 1.6635061e+01 1.23e-03 1.57e-02 -5.1 5.84e-02 - 1.00e+00 5.00e-01h 2
26r 1.6625746e+01 9.33e-04 1.18e-02 -5.1 4.01e-02 - 1.00e+00 2.50e-01h 3
27r 1.6596508e+01 9.06e-05 2.82e-07 -5.1 3.17e-02 - 1.00e+00 1.00e+00h 1
28r 1.6596381e+01 9.22e-09 8.47e-11 -5.1 2.61e-04 - 1.00e+00 1.00e+00h 1
29 1.6594177e+01 1.14e-05 2.02e+00 -8.6 4.33e+06 - 1.84e-08 1.23e-09f 4
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 1.6592864e+01 1.55e-05 2.02e+00 -8.6 5.04e+06 - 2.09e-07 1.19e-09f 5
31 1.6591729e+01 1.55e-05 2.02e+00 -8.6 1.40e+08 - 4.32e-08 1.42e-10f 4
32 1.6591237e+01 1.55e-05 2.02e+00 -8.6 2.05e+08 - 3.80e-08 4.25e-11f 5
33 1.6591121e+01 1.55e-05 2.02e+00 -8.6 1.74e+08 - 2.08e-08 1.17e-11f 7
34 1.6583523e+01 4.68e-04 2.02e+00 -8.6 2.95e+07 - 4.62e-07 4.39e-09f 1
35 1.6565886e+01 7.88e-04 2.02e+00 -8.6 4.76e+05 - 3.44e-05 6.52e-07f 1
36 1.6564668e+01 1.10e-01 2.02e+00 -8.6 2.61e+05 - 5.23e-05 1.91e-05f 1
37 1.6564309e+01 1.10e-01 2.02e+00 -8.6 1.97e+04 - 8.40e-04 4.89e-05h 1
38 1.6563516e+01 1.09e-01 2.02e+00 -8.6 4.32e+03 - 6.30e-03 9.65e-04h 1
39 1.6562582e+01 1.09e-01 2.01e+00 -8.6 2.63e+01 - 1.11e-02 7.61e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 1.6562162e+01 1.97e+00 2.00e+00 -8.6 2.11e+03 - 1.07e-03 3.31e-03h 1
41 1.6559812e+01 1.93e+00 1.96e+00 -8.6 3.28e+00 - 7.51e-05 2.05e-02h 1
42 1.6558871e+01 4.03e+00 1.94e+00 -8.6 5.95e+02 - 1.64e-03 6.37e-03h 1
In iteration 42, 1 Slack too small, adjusting variable bound
43 1.6558564e+01 4.03e+00 2.00e+00 -8.6 1.83e+02 - 3.83e-01 3.26e-03h 1
44 1.6527344e+01 2.71e+00 1.39e+00 -8.6 5.69e+00 - 9.10e-01 3.31e-01h 1
45 1.6473970e+01 1.64e-01 8.39e-02 -8.6 2.78e+00 - 1.00e+00 9.39e-01h 1
46 1.6470403e+01 1.69e-05 2.38e-03 -8.6 1.75e-01 - 9.83e-01 1.00e+00h 1
47 1.6470356e+01 5.52e-07 8.74e-06 -8.6 1.23e-02 -3.1 1.00e+00 1.00e+00h 1
48 1.6470356e+01 7.25e-09 1.84e-07 -8.6 1.18e-03 - 1.00e+00 1.00e+00h 1
In iteration 48, 1 Slack too small, adjusting variable bound
49 1.6470356e+01 4.13e-01 2.01e-01 -8.6 2.36e+03 - 1.03e-03 3.72e-03f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 1.6470356e+01 2.87e-01 9.98e-03 -8.6 8.42e-02 -3.6 1.00e+00 3.06e-01h 1
51 1.6470356e+01 1.43e-01 1.44e-03 -8.6 3.53e-01 - 9.93e-01 5.00e-01h 2
52 1.6470326e+01 1.47e-01 1.83e-01 -8.6 5.85e+02 - 7.78e-04 5.65e-03h 1
53 1.6470327e+01 1.42e-01 1.21e-02 -8.6 3.04e-02 -4.1 1.00e+00 3.56e-02h 1
54 1.6470355e+01 1.16e-04 7.92e-06 -8.6 2.93e-02 -4.6 1.00e+00 1.00e+00h 1
55 1.6470355e+01 1.10e-04 1.17e-02 -8.6 3.11e-01 - 9.37e-01 1.25e-01h 4
56 1.6470355e+01 1.07e-04 4.13e-02 -8.6 3.74e+01 - 7.41e-02 7.60e-04h 9
57 1.6470355e+01 1.01e-04 3.50e-02 -8.6 2.18e+00 - 5.37e-01 1.56e-02h 7
58 1.6470355e+01 1.00e-04 7.33e-02 -8.6 2.98e+01 - 9.67e-02 3.15e-04h 11
59 1.6470355e+01 9.20e-05 6.36e-02 -8.6 2.59e+00 - 2.40e-01 1.56e-02h 7
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 1.6470356e+01 2.32e-06 1.86e-07 -8.6 5.97e-03 -5.1 1.00e+00 1.00e+00h 1
61 1.6470356e+01 3.80e-05 4.28e-02 -8.6 1.22e+00 - 8.80e-01 6.25e-02h 5
62 1.6470356e+01 5.84e-05 8.40e-02 -8.6 2.64e+00 - 6.29e-01 3.12e-02h 6
63 1.6470356e+01 4.72e-06 3.78e-07 -8.6 1.48e-02 -5.5 1.00e+00 1.00e+00h 1
64 1.6470356e+01 3.78e-05 4.88e-02 -8.6 1.23e+00 - 9.68e-01 6.25e-02h 5
65 1.6470356e+01 2.74e-05 1.16e-02 -8.6 4.77e-02 -6.0 1.00e+00 2.75e-01h 2
66 1.6470356e+01 6.37e-05 4.89e-02 -8.6 1.23e+00 - 9.64e-01 6.25e-02h 5
67 1.6470356e+01 6.63e-05 4.73e-02 -8.6 3.26e+02 - 5.99e-04 1.83e-04h 8
68 1.6470356e+01 8.03e-05 6.43e-06 -8.6 6.10e-02 -6.5 1.00e+00 1.00e+00h 1
69 1.6470356e+01 7.76e-05 4.88e-02 -8.6 1.16e+00 - 9.52e-01 3.12e-02h 6
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 1.6470356e+01 7.57e-05 8.33e-02 -8.6 8.53e+01 - 6.21e-02 7.12e-04h 8
71 1.6470356e+01 8.90e-05 1.19e-01 -8.6 3.21e+01 - 2.73e-02 9.19e-04h 9
72 1.6470356e+01 9.53e-05 1.26e-01 -8.6 1.75e+02 - 3.06e-03 1.44e-04h 9
73 1.6470356e+01 9.66e-05 1.45e-01 -8.6 1.95e+02 - 1.33e-03 4.82e-05h 11
74 1.6470356e+01 2.62e-04 5.89e-05 -8.6 1.50e-01 -7.0 1.00e+00 1.00e+00H 1
75 1.6470356e+01 6.83e-03 3.75e-03 -8.6 1.00e+00 - 9.14e-01 1.00e+00h 1
76 1.6470356e+01 1.68e-04 1.34e-05 -8.6 8.54e-02 -7.4 1.00e+00 1.00e+00h 1
77 1.6470356e+01 3.25e-04 2.60e-05 -8.6 1.19e-01 -7.9 1.00e+00 1.00e+00h 1
78 1.6470356e+01 1.01e-02 9.28e-03 -8.6 1.59e+00 - 8.12e-01 1.00e+00h 1
79 1.6470356e+01 3.01e-04 2.39e-05 -8.6 1.20e-01 -8.4 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.6470356e+01 2.14e-04 5.85e-03 -8.6 3.17e-01 -8.9 1.00e+00 5.00e-01h 2
81 1.6470356e+01 6.56e-05 5.73e-06 -8.6 1.21e-01 -8.4 1.00e+00 1.00e+00h 1
82 1.6470356e+01 2.50e-04 4.38e-02 -8.6 3.59e+00 - 8.58e-01 6.25e-02h 5
In iteration 82, 1 Slack too small, adjusting variable bound
83 1.6470356e+01 6.39e-01 3.71e-02 -8.6 1.37e+01 - 1.00e+00 9.00e-01h 1
84 1.6470337e+01 2.99e-03 1.91e-04 -8.6 1.51e-01 -7.1 1.00e+00 1.00e+00h 1
85 1.6470356e+01 8.23e-06 7.09e-07 -8.6 4.07e-02 -7.6 1.00e+00 1.00e+00h 1
86 1.6470356e+01 1.66e-04 6.02e-02 -8.6 2.68e+00 - 7.94e-01 6.25e-02h 5
87 1.6470356e+01 1.66e-04 1.48e-01 -8.6 4.52e+00 - 1.00e+00 9.77e-04h 11
88 1.6470356e+01 1.66e-04 9.68e-02 -8.6 2.87e+00 - 8.30e-01 1.95e-03h 10
89 1.6470356e+01 1.66e-04 1.47e-01 -8.6 4.50e+00 - 1.00e+00 9.77e-04h 11
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 1.6470356e+01 1.66e-04 9.62e-02 -8.6 2.87e+00 - 8.34e-01 1.95e-03h 10
91 1.6470356e+01 1.66e-04 1.46e-01 -8.6 4.49e+00 - 1.00e+00 9.77e-04h 11
92 1.6470356e+01 1.65e-04 9.59e-02 -8.6 2.86e+00 - 8.36e-01 1.95e-03h 10
93 1.6470356e+01 1.65e-04 1.46e-01 -8.6 4.48e+00 - 1.00e+00 9.77e-04h 11
94 1.6470356e+01 1.65e-04 9.57e-02 -8.6 2.87e+00 - 8.39e-01 1.95e-03h 10
95 1.6470356e+01 1.65e-04 1.46e-01 -8.6 4.47e+00 - 1.00e+00 9.77e-04h 11
96 1.6470356e+01 4.68e-02 1.35e-02 -8.6 2.87e+00 - 8.42e-01 1.00e+00w 1
97 1.6470356e+01 1.45e-05 8.58e-07 -8.6 4.48e-02 -8.1 1.00e+00 1.00e+00h 1
98 1.6470356e+01 1.24e-05 5.24e-02 -8.6 1.70e+00 - 1.00e+00 1.56e-02h 7
99 1.6470356e+01 1.26e-05 5.71e-02 -8.6 1.81e+00 - 1.00e+00 7.81e-03h 8
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 1.6470356e+01 1.27e-05 5.41e-02 -8.6 1.72e+00 - 1.00e+00 3.91e-03h 9
101 1.6470356e+01 1.27e-05 5.49e-02 -8.6 1.74e+00 - 1.00e+00 3.91e-03h 9
102 1.6470356e+01 1.27e-05 5.43e-02 -8.6 1.73e+00 - 1.00e+00 3.91e-03h 9
103 1.6470356e+01 1.27e-05 5.45e-02 -8.6 1.73e+00 - 1.00e+00 1.95e-03h 10
104 1.6470356e+01 1.27e-05 5.44e-02 -8.6 1.72e+00 - 1.00e+00 9.77e-04h 11
105 1.6470356e+01 1.27e-05 5.44e-02 -8.6 1.72e+00 - 1.00e+00 9.77e-04h 11
106 1.6470356e+01 1.27e-05 5.43e-02 -8.6 1.72e+00 - 1.00e+00 9.77e-04h 11
107 1.6470356e+01 1.27e-05 5.43e-02 -8.6 1.72e+00 - 1.00e+00 9.77e-04h 11
108 1.6470356e+01 1.57e-02 1.22e-03 -8.6 1.72e+00 - 1.00e+00 1.00e+00w 1
109 1.6470356e+01 8.19e-07 6.54e-08 -8.6 1.84e-02 -8.5 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 1.6470356e+01 7.52e-07 1.00e-02 -8.6 3.64e-01 - 1.00e+00 3.12e-02h 6
111 1.6470356e+01 7.78e-07 9.85e-03 -8.6 3.52e-01 - 1.00e+00 1.56e-02h 7
112 1.6470356e+01 7.78e-07 9.80e-03 -8.6 3.46e-01 - 1.00e+00 3.91e-03h 9
113 1.6470356e+01 7.78e-07 9.77e-03 -8.6 3.45e-01 - 1.00e+00 3.91e-03h 9
114 1.6470356e+01 7.77e-07 9.75e-03 -8.6 3.44e-01 - 1.00e+00 1.95e-03h 10
115 1.6470356e+01 7.76e-07 9.73e-03 -8.6 3.43e-01 - 1.00e+00 1.95e-03h 10
116 1.6470356e+01 7.75e-07 9.71e-03 -8.6 3.42e-01 - 1.00e+00 1.95e-03h 10
117 1.6470356e+01 7.75e-07 9.69e-03 -8.6 3.42e-01 - 1.00e+00 1.95e-03h 10
118 1.6470356e+01 7.74e-07 9.68e-03 -8.6 3.41e-01 - 1.00e+00 1.95e-03h 10
119 1.6470356e+01 7.73e-07 9.66e-03 -8.6 3.40e-01 - 1.00e+00 1.95e-03h 10
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
120 1.6470356e+01 5.77e-04 4.92e-05 -8.6 3.40e-01 - 1.00e+00 1.00e+00w 1
121 1.6470356e+01 8.19e-10 5.80e-11 -8.6 2.13e-04 -8.1 1.00e+00 1.00e+00h 1
In iteration 121, 1 Slack too small, adjusting variable bound
122 1.6470356e+01 2.30e-02 4.10e-02 -10.0 3.52e+00 - 5.96e-01 7.65e-01h 1
123 1.6470351e+01 8.72e-03 7.11e-04 -10.0 5.18e-01 -8.6 1.00e+00 1.00e+00h 1
124 1.6470353e+01 8.75e-03 7.82e-02 -10.0 6.45e-01 -9.1 1.00e+00 5.00e-01h 2
125 1.6470354e+01 4.45e-03 1.01e-02 -10.0 7.00e-02 -8.7 1.00e+00 5.00e-01h 2
126 1.6470355e+01 2.24e-03 6.27e-03 -10.0 1.59e-01 -9.1 1.00e+00 5.00e-01h 2
127 1.6470355e+01 1.14e-03 7.54e-03 -10.0 3.56e-01 -9.6 1.00e+00 5.00e-01h 2
128 1.6470355e+01 1.14e-03 4.27e-01 -10.0 2.44e+01 - 3.49e-01 2.32e-04h 12
129 1.6470355e+01 1.14e-03 1.20e+00 -10.0 2.42e+01 - 1.00e+00 2.58e-04h 12
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
130 1.6470355e+01 1.14e-03 1.04e+00 -10.0 1.84e+01 - 4.85e-01 6.79e-04h 11
131 1.6470355e+01 1.14e-03 8.52e-01 -10.0 1.81e+01 - 1.00e+00 6.88e-04h 11
132 1.6470355e+01 1.14e-03 7.67e-01 -10.0 1.56e+01 - 5.96e-01 8.07e-04h 11
133 1.6470355e+01 1.14e-03 7.20e-01 -10.0 1.57e+01 - 1.00e+00 8.03e-04h 11
134 1.6470356e+01 6.63e-01 6.52e-02 -10.0 1.45e+01 - 6.48e-01 9.05e-01w 1
135 1.6469921e+01 4.99e-01 4.85e-01 -10.0 2.31e+00 -6.5 9.39e-01 3.23e-01h 1
136 1.6470010e+01 2.84e-01 1.71e-01 -10.0 4.52e-01 -6.0 6.31e-01 4.66e-01h 1
137 1.6470271e+01 5.04e-02 2.74e-02 -10.0 4.55e-01 -6.5 1.00e+00 8.22e-01h 1
138 1.6470355e+01 5.07e-04 4.07e-05 -10.0 6.62e-02 -7.0 1.00e+00 1.00e+00h 1
139 1.6470355e+01 5.07e-04 8.10e-02 -10.0 1.99e+00 - 2.37e-01 3.91e-03h 9
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
140 1.6470355e+01 5.06e-04 3.09e-01 -10.0 2.86e+00 - 6.66e-01 1.95e-03h 10
141 1.6470355e+01 5.06e-04 3.15e-01 -10.0 4.43e+00 - 1.80e-01 9.77e-04h 11
142 1.6470355e+01 5.06e-04 3.00e-01 -10.0 4.12e+00 - 1.25e-01 9.77e-04h 11
143 1.6470355e+01 5.06e-04 2.98e-01 -10.0 5.45e+00 - 6.05e-01 4.88e-04h 12
144 1.6470355e+01 5.06e-04 5.09e-01 -10.0 1.60e+01 - 2.86e-02 1.08e-04h 6
145 1.6470355e+01 5.06e-04 4.07e-01 -10.0 3.54e+00 - 2.07e-01 9.77e-04h 11
146 1.6470355e+01 5.05e-04 5.06e-01 -10.0 5.62e+00 - 1.00e+00 4.88e-04h 12
147 1.6470355e+01 5.05e-04 3.92e-01 -10.0 3.59e+00 - 1.51e-01 9.40e-04h 9
148 1.6470355e+01 5.05e-04 3.39e-01 -10.0 3.46e+00 - 2.11e-01 9.77e-04h 11
149 1.6470356e+01 3.87e-01 1.55e-01 -10.0 5.18e+00 - 2.49e-01 1.00e+00w 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
150 1.6470217e+01 3.12e-01 6.78e-02 -10.0 4.36e+00 -7.5 5.41e-01 5.20e-01h 1
151 1.6470254e+01 2.34e-01 8.96e-02 -10.0 7.76e-01 -7.1 1.00e+00 2.50e-01h 3
152 1.6470269e+01 1.96e-01 4.38e-02 -10.0 3.09e-01 -7.5 1.00e+00 1.63e-01h 3
153 1.6470351e+01 1.20e-03 9.41e-05 -10.0 5.77e-02 -7.1 1.00e+00 1.00e+00h 1
154 1.6470356e+01 1.35e-04 2.26e-03 -10.0 3.59e-02 -7.6 8.29e-01 1.00e+00h 1
155 1.6470356e+01 1.70e-04 1.38e-01 -10.0 3.40e+00 - 2.29e-01 1.56e-02h 7
156 1.6470356e+01 3.14e-04 1.69e-01 -10.0 5.24e+00 - 1.60e-01 1.26e-02h 7
157 1.6470356e+01 3.13e-04 2.03e-01 -10.0 4.78e+00 - 1.39e-01 2.44e-04h 13
158 1.6470356e+01 3.13e-04 4.60e-01 -10.0 5.82e+00 - 3.21e-01 2.44e-04h 13
159 1.6470356e+01 3.13e-04 4.39e-01 -10.0 4.70e+00 - 1.54e-01 4.88e-04h 12
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
160 1.6470356e+01 3.13e-04 5.90e-01 -10.0 5.87e+00 - 2.72e-01 2.44e-04h 13
161 1.6470356e+01 3.13e-04 5.57e-01 -10.0 4.80e+00 - 1.47e-01 2.44e-04h 13
162 1.6470356e+01 3.13e-04 6.86e-01 -10.0 5.87e+00 - 2.97e-01 2.44e-04h 13
163 1.6470356e+01 3.13e-04 6.35e-01 -10.0 4.77e+00 - 1.51e-01 2.44e-04h 13
164 1.6470356e+01 3.13e-04 7.35e-01 -10.0 5.87e+00 - 2.87e-01 2.44e-04h 13
165 1.6470356e+01 4.94e-01 3.59e-01 -10.0 4.80e+00 - 1.50e-01 1.00e+00w 1
166 1.6470188e+01 3.81e-01 7.49e-02 -10.0 2.59e+00 -7.2 9.26e-01 7.16e-01h 1
167 1.6470257e+01 1.65e-01 2.14e-01 -10.0 7.16e-01 -6.7 1.00e+00 5.00e-01h 2
168 1.6470346e+01 6.51e-03 5.38e-04 -10.0 1.18e-01 -7.2 1.00e+00 1.00e+00h 1
169 1.6470356e+01 3.88e-04 6.98e-04 -10.0 7.98e-02 -7.7 9.74e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
170 1.6470356e+01 3.88e-04 1.47e-01 -10.0 2.35e+00 - 3.34e-01 9.77e-04h 11
171 1.6470356e+01 3.88e-04 4.32e-01 -10.0 2.38e+00 - 1.00e+00 9.77e-04h 11
172 1.6470356e+01 3.87e-04 4.29e-01 -10.0 2.32e+00 - 3.36e-01 9.77e-04h 11
173 1.6470356e+01 3.87e-04 4.30e-01 -10.0 2.39e+00 - 1.00e+00 9.77e-04h 11
174 1.6470356e+01 3.87e-04 4.27e-01 -10.0 2.32e+00 - 3.38e-01 9.77e-04h 11
175 1.6470356e+01 3.87e-04 4.28e-01 -10.0 2.39e+00 - 1.00e+00 9.77e-04h 11
176 1.6470356e+01 3.87e-04 4.26e-01 -10.0 2.32e+00 - 3.40e-01 9.77e-04h 11
177 1.6470356e+01 3.86e-04 4.27e-01 -10.0 2.38e+00 - 1.00e+00 9.77e-04h 11
178 1.6470356e+01 3.86e-04 4.25e-01 -10.0 2.32e+00 - 3.42e-01 9.77e-04h 11
179 1.6470356e+01 3.86e-04 4.26e-01 -10.0 2.38e+00 - 1.00e+00 9.77e-04h 11
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
180 1.6470356e+01 2.59e-01 2.46e-02 -10.0 2.32e+00 - 3.43e-01 1.00e+00w 1
181 1.6470349e+01 6.18e-03 5.06e-04 -10.0 1.12e-01 -7.3 1.00e+00 1.00e+00h 1
182 1.6470356e+01 1.71e-05 1.36e-06 -10.0 2.59e-02 -7.7 1.00e+00 1.00e+00h 1
183 1.6470356e+01 1.10e-04 1.41e-01 -10.0 8.49e-01 - 1.00e+00 6.25e-02h 5
184 1.6470356e+01 1.96e-04 1.42e-01 -10.0 8.71e-01 - 1.00e+00 6.25e-02h 5
185 1.6470356e+01 2.12e-04 1.45e-01 -10.0 8.81e-01 - 1.00e+00 3.12e-02h 6
186 1.6470356e+01 2.27e-04 1.41e-01 -10.0 8.69e-01 - 1.00e+00 3.12e-02h 6
187 1.6470356e+01 2.41e-04 1.38e-01 -10.0 8.55e-01 - 1.00e+00 3.12e-02h 6
188 1.6470356e+01 2.02e-02 1.61e-03 -10.0 8.40e-01 - 1.00e+00 1.00e+00h 1
189 1.6470356e+01 1.24e-04 1.02e-05 -10.0 1.73e-02 -8.2 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
190 1.6470356e+01 3.26e-01 1.01e-01 -10.0 3.15e-01 - 2.18e-02 3.21e-02H 1
191 1.6470356e+01 3.26e-01 1.81e-01 -10.0 2.12e+00 -6.9 1.00e+00 2.44e-04h 13
192 1.6470306e+01 1.37e-02 4.77e-03 -10.0 4.62e-01 -7.4 9.76e-01 1.00e+00h 1
193 1.6470331e+01 6.37e-03 7.40e-02 -10.0 1.83e+00 -7.8 1.00e+00 5.00e-01h 2
194 1.6470343e+01 3.19e-03 5.62e-03 -10.0 1.37e-01 -8.3 9.99e-01 5.00e-01h 2
195 1.6470350e+01 1.60e-03 1.44e-02 -10.0 3.82e-01 -8.8 1.00e+00 5.00e-01h 2
196 1.6470350e+01 1.40e-03 6.14e-02 -10.0 1.00e+00 -9.3 1.00e+00 1.25e-01h 4
197 1.6470350e+01 1.39e-03 1.90e-01 -10.0 3.29e+00 -9.7 1.00e+00 1.56e-02h 7
198 1.6470351e+01 1.23e-03 6.32e-02 -10.0 1.06e+00 -9.3 1.00e+00 1.25e-01h 4
199 1.6470351e+01 1.23e-03 1.70e-01 -10.0 3.41e+00 -9.8 8.47e-01 1.56e-02h 7
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
200 1.6470352e+01 1.26e-03 6.47e-02 -10.0 1.12e+00 -9.4 1.00e+00 1.25e-01h 4
201 1.6470352e+01 1.26e-03 1.61e-01 -10.0 3.56e+00 -9.8 7.60e-01 1.56e-02h 7
202 1.6470352e+01 1.30e-03 6.60e-02 -10.0 1.17e+00 -9.4 1.00e+00 1.25e-01h 4
203 1.6470356e+01 1.39e-01 5.06e-02 -10.0 3.78e+00 -9.9 7.04e-01 1.00e+00w 1
204 1.6470347e+01 2.63e-03 2.07e-04 -10.0 2.12e-01 -7.7 1.00e+00 1.00e+00h 1
205 1.6470356e+01 2.98e-05 2.37e-06 -10.0 4.40e-02 -8.1 1.00e+00 1.00e+00h 1
206 1.6470356e+01 7.48e-05 1.26e-01 -10.0 2.64e+00 - 1.00e+00 3.12e-02h 6
207 1.6470356e+01 8.77e-05 1.26e-01 -10.0 2.60e+00 - 1.00e+00 1.56e-02h 7
208 1.6470356e+01 9.04e-05 1.24e-01 -10.0 2.54e+00 - 1.00e+00 7.81e-03h 8
209 1.6470356e+01 9.30e-05 1.23e-01 -10.0 2.52e+00 - 1.00e+00 7.81e-03h 8
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
210 1.6470356e+01 9.35e-05 1.22e-01 -10.0 2.49e+00 - 1.00e+00 3.91e-03h 9
211 1.6470356e+01 5.27e-02 4.38e-03 -10.0 2.48e+00 - 1.00e+00 1.00e+00h 1
212 1.6470356e+01 2.91e-05 2.16e-06 -10.0 2.08e-02 -8.6 1.00e+00 1.00e+00h 1
213 1.6470356e+01 2.90e-05 2.84e-02 -10.0 6.07e-01 - 1.00e+00 7.81e-03h 8
214 1.6470356e+01 2.90e-05 2.82e-02 -10.0 6.03e-01 - 1.00e+00 7.81e-03h 8
215 1.6470356e+01 2.89e-05 2.80e-02 -10.0 5.98e-01 - 1.00e+00 7.81e-03h 8
216 1.6470356e+01 2.89e-05 2.78e-02 -10.0 5.93e-01 - 1.00e+00 7.81e-03h 8
217 1.6470356e+01 2.88e-05 2.76e-02 -10.0 5.88e-01 - 1.00e+00 7.81e-03h 8
218 1.6470356e+01 2.88e-05 2.74e-02 -10.0 5.84e-01 - 1.00e+00 7.81e-03h 8
219 1.6470356e+01 2.87e-05 2.71e-02 -10.0 5.79e-01 - 1.00e+00 7.81e-03h 8
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
220 1.6470356e+01 2.87e-05 2.69e-02 -10.0 5.75e-01 - 1.00e+00 7.81e-03h 8
221 1.6470356e+01 2.86e-05 2.67e-02 -10.0 5.70e-01 - 1.00e+00 7.81e-03h 8
222 1.6470356e+01 2.86e-05 2.65e-02 -10.0 5.66e-01 - 1.00e+00 7.81e-03h 8
223 1.6470356e+01 2.63e-03 2.18e-04 -10.0 5.61e-01 - 1.00e+00 1.00e+00w 1
224 1.6470356e+01 3.26e-07 2.45e-08 -10.0 2.15e-03 -9.1 1.00e+00 1.00e+00h 1
225 1.6470356e+01 3.07e-07 5.29e-04 -10.0 2.43e-02 - 1.00e+00 6.25e-02h 5
226 1.6470356e+01 2.88e-07 4.96e-04 -10.0 2.28e-02 - 1.00e+00 6.25e-02h 5
227 1.6470356e+01 2.70e-07 4.65e-04 -10.0 2.14e-02 - 1.00e+00 6.25e-02h 5
228 1.6470356e+01 2.38e-07 4.06e-04 -10.0 2.00e-02 - 1.00e+00 1.25e-01h 4
229 1.6470356e+01 2.10e-07 3.56e-04 -10.0 1.75e-02 - 1.00e+00 1.25e-01h 4
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
230 1.6470356e+01 1.85e-07 3.11e-04 -10.0 1.53e-02 - 1.00e+00 1.25e-01h 4
231 1.6470356e+01 1.41e-07 2.33e-04 -10.0 1.34e-02 - 1.00e+00 2.50e-01h 3
232 1.6470356e+01 1.15e-07 1.75e-04 -10.0 1.01e-02 - 1.00e+00 2.50e-01h 3
233 1.6470356e+01 1.21e-07 8.74e-05 -10.0 7.54e-03 - 1.00e+00 5.00e-01h 2
234 1.6470356e+01 6.29e-08 5.25e-09 -10.0 3.77e-03 - 1.00e+00 1.00e+00h 1
235 1.6470356e+01 1.42e-14 1.07e-15 -10.0 1.24e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 235
(scaled) (unscaled)
Objective...............: 1.6470355768670274e+01 1.6470355768670274e+01
Dual infeasibility......: 1.0698114603874582e-15 1.0698114603874582e-15
Constraint violation....: 1.4210854715202004e-14 1.4210854715202004e-14
Complementarity.........: 9.0909133698307485e-11 9.0909133698307485e-11
Overall NLP error.......: 9.0909133698307485e-11 9.0909133698307485e-11
Number of objective function evaluations = 1395
Number of objective gradient evaluations = 221
Number of equality constraint evaluations = 1395
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 237
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 235
Total CPU secs in IPOPT (w/o function evaluations) = 0.063
Total CPU secs in NLP function evaluations = 0.010
EXIT: Optimal Solution Found.
Moving elements: [1.771518472796034, 1.7709493825300293, 1.783521573771728, 1.7956856441920985, 1.8050747071322701, 1.8115827813436336, 1.8154385505809003, 0.9076421104172527, 0.9076421567549828, 2.101300389151346]
Controls: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -3.0, -3.0, -3.0]
Final time: 16.470355768670274
Check the switch and continuous trajectory¶
For this model the bang-bang reference accelerates at until speed , then brakes at . Integrating gives
Solve with . The two durations are
These formulas come from the continuous ODE, independently of collocation. They define a reference candidate, not a general global-optimality certificate for arbitrary nonlinear control problems.
The next cell also integrates each computed piecewise-constant control on its actual intervals using a tight-tolerance ODE solve, restarting exactly at each control change. This checks the continuous terminal residual separately from the collocation constraints.
On a uniform mesh, a boundary where the control first changes sign is not a true bang-bang switch if an intermediate-control element remains. Read that interval and its control value. The moving-element run must put controls at the acceleration/braking bounds on either side before its boundary can be compared with the one-switch reference. Several length allocations within a constant-control arc can be nearly equivalent, so do not interpret the reported mesh as unique.
from scipy.optimize import brentq
from scipy.integrate import solve_ivp
import scipy, pyomo, platform
R, acceleration, braking, distance = 0.001, 1.0, 3.0, 100.0
v_switch = brentq(lambda v: (-np.log1p(-R*v*v/acceleration)
+np.log1p(R*v*v/braking))/(2*R)-distance,
0, np.sqrt(acceleration/R)*(1-1e-12))
t_switch = np.arctanh(v_switch*np.sqrt(R/acceleration))/np.sqrt(acceleration*R)
t_brake = np.arctan(v_switch*np.sqrt(R/braking))/np.sqrt(braking*R)
reference = {"switch_speed":v_switch,"switch_time":t_switch,"final_time":t_switch+t_brake}
uniform = [extract_race_car_results(m_coarse,NFE_COARSE),
extract_race_car_results(m_fine,NFE_FINE)]
edges = np.r_[0.0,np.cumsum([pyo.value(m_variable.h[i]) for i in m_variable.I])]
variable = {"label":"N = 10, moving","nfe":10,"tf":pyo.value(m_variable.tf),
"elements":edges.tolist(),
"element_u":[pyo.value(m_variable.u[i]) for i in m_variable.I]}
comparisons = []
for mesh in uniform+[variable]:
state = np.array([0.0,0.0])
for left,right,control in zip(mesh["elements"][:-1],mesh["elements"][1:],mesh["element_u"]):
simulation = solve_ivp(lambda t,y:[y[1],control-R*y[1]**2],
(left,right),state,rtol=1e-11,atol=1e-12,method="DOP853")
assert simulation.success and np.isclose(simulation.t[-1],right)
state = simulation.y[:,-1]
switches = [mesh["elements"][j] for j in range(1,len(mesh["element_u"]))
if mesh["element_u"][j-1]>0 and mesh["element_u"][j]<0]
comparisons.append({**mesh,"simulated_x_final":float(state[0]),
"simulated_v_final":float(state[1]),
"switch_boundaries":switches})
moving_check = comparisons[-1]
assert all(min(abs(u-acceleration),abs(u+braking))<1e-6 for u in variable["element_u"])
assert len(moving_check["switch_boundaries"])==1
assert abs(moving_check["switch_boundaries"][0]-reference["switch_time"])<1e-5
assert abs(moving_check["tf"]-reference["final_time"])<1e-5
assert abs(moving_check["simulated_x_final"]-distance)<1e-5
assert abs(moving_check["simulated_v_final"])<1e-6
comparison_results = {"meshes":comparisons,"reference":reference,"R":R,"L":distance,
"u_lower":-3.0,"u_upper":1.0,
"environment":{"pyomo":pyomo.__version__,"scipy":scipy.__version__,
"python":platform.python_version()}}
for row in comparisons:
print(f'{row["label"]:>16}: tf={row["tf"]:.8f}, '
f'continuous dx={row["simulated_x_final"]-distance:+.3e}, '
f'continuous v={row["simulated_v_final"]:+.3e}, '
f'switch boundaries={row["switch_boundaries"]}')
print("Independent reference:",reference)
helper.save_results("race-car-variable-elements",comparison_results,
notebook="notebooks/3-dev/PyomoDAE_car.ipynb",
source_tag="handout:race-car-model",
source_cells=[("notebooks/3-dev/PyomoDAE_car.ipynb",tag) for tag in
("handout-output:race-car-coarse","handout-output:race-car-fine",
"numeric:race-car-extract","numeric:race-car-variable",
"numeric:race-car-independent")],
description="Uniform10/100 and optimized10-element car controls, with independent bang-bang and continuous-integration checks.",
solver="Ipopt3.13.2/ma27; DOP853 independent trajectory check")
results = comparison_results
N = 10: tf=16.57829696, continuous dx=-1.367e-06, continuous v=+3.260e-08, switch boundaries=[11.604807872512213]
N = 100: tf=16.47113797, continuous dx=+3.752e-12, continuous v=+1.502e-12, switch boundaries=[12.518064854519096]
N = 10, moving: tf=16.47035577, continuous dx=-3.350e-06, continuous v=+4.167e-08, switch boundaries=[12.553771112346693]
Independent reference: {'switch_speed': 11.933371838610697, 'switch_time': np.float64(12.55377133285438), 'final_time': np.float64(16.470356063247646)}
[helper] wrote figures/results/race-car-variable-elements.json
def plot_variable_element_car(results):
fig, axes = plt.subplots(1,2,figsize=(10.5,4.2),gridspec_kw={"width_ratios":[1.25,1]})
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
styles = ["-","--","-."]
for j,row in enumerate(results["meshes"]):
axes[0].stairs(row["element_u"],row["elements"],baseline=None,
label=row["label"],linestyle=styles[j],color=colors[j],linewidth=2)
if j==2:
for edge in row["elements"][1:-1]:
axes[0].axvline(edge,color="0.8",linewidth=.6,zorder=0)
axes[0].axvline(results["reference"]["switch_time"],color="black",
linestyle=":",label="continuous switch",linewidth=1.5)
axes[0].set(xlabel="Physical time (s)",ylabel="Control (m/s²)",ylim=(-3.3,1.35))
axes[0].legend(loc="lower left",fontsize=10)
moving = results["meshes"][-1]
lengths = np.diff(moving["elements"])
axes[1].bar(np.arange(1,len(lengths)+1),lengths,color=colors[2],
edgecolor="black",hatch="//",alpha=.65)
axes[1].set(xlabel="Finite element",ylabel="Optimized length (s)")
axes[1].set_xticks(np.arange(1,len(lengths)+1))
fig.tight_layout()
return fig
fig = plot_variable_element_car(results)
helper.save_figure(fig,"race-car-variable-elements")
[helper] wrote media/figures/race-car-variable-elements.png and .pdf
# Same physical duration bound in every row; controls are not reduced here.
from tempfile import TemporaryDirectory
from pathlib import Path
import re
initialization_results = []
for requested_tf in [1.0,0.0,1e6]:
trial = create_model1()
for state_var in [trial.x,trial.v]:
for v in state_var.values():v.set_value(0.0)
trial.tf.set_value(requested_tf,skip_validation=True)
pyo.TransformationFactory("dae.collocation").apply_to(trial,nfe=15,ncp=3,
scheme="LAGRANGE-RADAU")
with TemporaryDirectory() as directory:
logfile = Path(directory)/"ipopt.log"
status = pyo.SolverFactory("ipopt").solve(trial, logfile=str(logfile))
transcript = logfile.read_text()
iteration_match = re.search(r"Number of Iterations.*?:\s*(\d+)", transcript)
assert iteration_match, "Missing iteration summary in actual Ipopt log"
assert pyo.check_optimal_termination(status)
record = {"initial_tf":requested_tf,"lower_tf":trial.tf.lb,
"final_tf":pyo.value(trial.tf),"iterations":int(iteration_match.group(1))}
initialization_results.append(record)
print(record)
print(transcript)
assert max(abs(row["final_tf"]-16.4708) for row in initialization_results)<1e-3
{'initial_tf': 1.0, 'lower_tf': 0.001, 'final_tf': 16.47083526594702, 'iterations': 81}
Solver command line: ['/Users/adowling/.idaes/bin/ipopt', '/var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmpz5kxj919.pyomo.nl', '-AMPL']
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 769
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 135
Total number of variables............................: 228
variables with only lower bounds: 47
variables with lower and upper bounds: 91
variables with only upper bounds: 0
Total number of equality constraints.................: 184
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 1.00e+02 4.11e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 4.3799862e+00 9.98e+01 1.18e+02 -1.0 2.29e+03 - 9.91e-05 1.47e-03f 1
2 4.6649766e+00 9.98e+01 1.18e+02 -1.0 5.92e+03 - 3.24e-04 3.99e-04h 1
3 4.7457534e+00 9.98e+01 1.18e+02 -1.0 2.89e+03 - 1.44e-03 1.76e-04h 1
4 4.9686989e+00 9.97e+01 1.17e+02 -1.0 2.08e+03 - 6.52e-04 5.10e-04h 1
5 5.5148707e+00 9.96e+01 1.16e+02 -1.0 1.82e+03 - 1.12e-03 1.30e-03h 1
6 5.6569588e+00 9.96e+01 1.16e+02 -1.0 1.60e+03 - 7.03e-04 3.76e-04h 1
7 6.0741393e+00 9.95e+01 1.15e+02 -1.0 1.58e+03 - 7.17e-04 1.13e-03h 2
8 6.4685847e+00 9.93e+01 1.14e+02 -1.0 1.55e+03 - 1.40e-03 1.13e-03h 1
9 6.9397559e+00 9.92e+01 1.13e+02 -1.0 1.43e+03 - 2.59e-03 1.42e-03h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 7.5785785e+00 9.90e+01 1.12e+02 -1.0 1.23e+03 - 9.95e-03 2.05e-03h 1
11 8.0307614e+00 9.88e+01 1.11e+02 -1.0 9.16e+02 - 7.86e-03 1.59e-03h 3
12 8.4973883e+00 9.87e+01 1.11e+02 -1.0 8.15e+02 - 3.86e-03 1.73e-03h 3
13 9.4357669e+00 9.83e+01 1.09e+02 -1.0 7.76e+02 - 7.23e-03 3.64e-03h 2
14 1.0354892e+01 9.79e+01 1.07e+02 -1.0 7.17e+02 - 5.04e-03 3.87e-03h 2
15 1.1142633e+01 9.76e+01 1.06e+02 -1.0 7.32e+02 - 5.40e-03 3.60e-03h 2
16 1.2248161e+01 9.70e+01 1.04e+02 -1.0 7.29e+02 - 5.20e-02 5.47e-03h 1
17 1.3492428e+01 9.64e+01 1.02e+02 -1.0 7.03e+02 - 2.06e-03 6.83e-03h 3
18 1.4386690e+01 9.59e+01 1.01e+02 -1.0 7.25e+02 - 1.34e-02 5.29e-03h 2
19 1.5212292e+01 9.54e+01 1.00e+02 -1.0 7.17e+02 - 1.47e-02 5.29e-03h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.6178110e+01 9.47e+01 1.00e+02 -1.0 7.13e+02 - 8.80e-02 6.59e-03h 2
21 1.7727818e+01 9.37e+01 1.28e+02 -1.0 7.15e+02 - 4.32e-02 1.09e-02h 3
22 1.9701570e+01 9.23e+01 1.23e+02 -1.0 7.20e+02 - 3.39e-03 1.46e-02h 3
23 2.1016777e+01 9.14e+01 1.23e+02 -1.0 7.38e+02 - 2.60e-02 1.03e-02h 3
24 2.2664055e+01 9.02e+01 1.21e+02 -1.0 7.42e+02 - 1.76e-02 1.33e-02h 2
25 2.4090998e+01 8.91e+01 1.16e+02 -1.0 7.49e+02 - 6.87e-02 1.20e-02h 1
26 2.5525143e+01 8.77e+01 1.44e+02 -1.0 6.82e+02 - 1.09e-01 1.51e-02h 1
27 2.7301499e+01 8.62e+01 1.35e+02 -1.0 7.00e+02 - 2.68e-03 1.82e-02h 2
28 2.7949095e+01 8.56e+01 1.53e+02 -1.0 7.32e+02 - 9.69e-02 6.48e-03h 1
29 2.9928433e+01 8.37e+01 1.73e+02 -1.0 6.78e+02 - 6.82e-02 2.26e-02h 3
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 3.2639039e+01 8.13e+01 1.70e+02 -1.0 7.40e+02 - 5.01e-02 2.86e-02h 3
31 3.6054331e+01 7.83e+01 1.81e+02 -1.0 7.77e+02 - 2.22e-01 3.66e-02h 2
32 3.8559981e+01 7.61e+01 1.73e+02 -1.0 8.11e+02 - 1.67e-02 2.77e-02h 2
33 4.0356631e+01 7.47e+01 1.69e+02 -1.0 8.82e+02 - 1.40e-01 1.90e-02h 2
34 4.2371939e+01 7.29e+01 1.69e+02 -1.0 8.27e+02 - 4.30e-02 2.40e-02h 1
35 4.4173719e+01 7.06e+01 1.64e+02 -1.0 6.38e+02 - 3.86e-02 3.19e-02h 1
36 4.7029410e+01 6.69e+01 1.54e+02 -1.0 6.39e+02 - 4.52e-02 5.19e-02h 3
37 5.1596967e+01 6.10e+01 1.60e+02 -1.0 6.31e+02 - 3.82e-02 8.84e-02h 1
38 5.3121881e+01 5.88e+01 1.25e+02 -1.0 5.73e+02 - 2.91e-01 3.60e-02h 1
39 5.5663687e+01 5.62e+01 1.20e+02 -1.0 7.34e+02 - 4.80e-02 4.33e-02h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 6.0647331e+01 5.07e+01 8.79e+01 -1.0 6.75e+02 - 3.42e-02 9.87e-02H 1
41 6.1107677e+01 4.99e+01 8.49e+01 -1.0 4.94e+02 - 2.91e-02 1.48e-02h 1
42 8.6175359e+01 1.77e+01 2.03e+02 -1.0 4.74e+02 - 5.38e-02 9.90e-01H 1
43 8.4641323e+01 3.34e-01 1.17e+03 -1.0 9.82e+00 - 4.06e-01 9.90e-01h 1
44 7.9245451e+01 5.51e+00 5.38e+04 -1.0 2.17e+02 - 1.25e-01 3.90e-01f 1
45 8.2367185e+01 2.09e+00 2.49e+05 -1.0 5.50e+01 - 2.30e-01 1.00e+00h 1
46 7.8890218e+01 3.45e+00 1.31e+06 -1.0 1.16e+02 - 5.91e-01 4.45e-01f 1
47 7.8281979e+01 5.02e-02 2.34e+06 -1.0 7.75e+00 - 4.22e-01 1.00e+00f 1
48 7.4953257e+01 3.21e+00 2.05e+06 -1.0 3.99e+02 - 1.24e-01 1.86e-01f 1
49 7.0169579e+01 5.75e+00 1.97e+06 -1.0 3.11e+02 - 3.61e-02 2.20e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 6.8227989e+01 4.15e+00 1.85e+06 -1.0 1.66e+02 - 6.15e-02 3.89e-01f 1
51 6.8348122e+01 3.82e+00 1.50e+06 -1.0 3.94e+02 - 1.89e-01 7.92e-02f 1
52 6.6615017e+01 3.19e+00 9.79e+05 -1.0 1.10e+02 - 3.49e-01 3.04e-01f 1
53 6.1180073e+01 8.70e+00 9.51e+05 -1.0 4.84e+02 - 2.84e-02 1.54e-01f 1
54 6.6622345e+01 4.84e+00 7.12e+05 -1.0 1.01e+02 - 2.52e-01 7.08e-01H 1
55 6.1498977e+01 8.10e+00 7.04e+05 -1.0 2.56e+02 - 1.11e-02 2.87e-01f 1
56 5.5676542e+01 1.10e+01 5.15e+05 -1.0 1.41e+02 - 2.69e-01 5.34e-01f 1
57 5.4189404e+01 1.07e+01 3.85e+05 -1.0 2.98e+02 - 2.52e-01 6.54e-02f 1
58 5.1089198e+01 9.38e+00 3.38e+05 -1.0 3.15e+02 - 1.22e-01 2.93e-01f 1
59 5.0086150e+01 8.12e+00 3.01e+05 -1.0 2.00e+02 - 1.09e-01 1.51e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 4.8374830e+01 7.69e+00 2.46e+05 -1.0 3.62e+02 - 1.85e-01 1.17e-01f 1
61 4.5650595e+01 6.74e+00 2.19e+05 -1.0 2.83e+02 - 1.07e-01 2.56e-01f 1
62 4.3274170e+01 6.37e+00 1.70e+05 -1.0 2.17e+02 - 2.24e-01 2.06e-01f 1
63 3.9655856e+01 7.31e+00 1.31e+05 -1.0 2.71e+02 - 2.28e-01 2.37e-01f 1
64 3.6828472e+01 7.48e+00 9.92e+04 -1.0 2.37e+02 - 2.45e-01 1.91e-01f 1
65 3.2021608e+01 1.10e+01 7.65e+04 -1.0 2.30e+02 - 2.28e-01 3.18e-01f 1
66 2.7379470e+01 1.20e+01 5.78e+04 -1.0 1.53e+02 - 2.44e-01 4.36e-01f 1
67 2.4026960e+01 1.06e+01 3.72e+04 -1.0 1.51e+02 - 3.57e-01 4.62e-01f 1
68 2.0195133e+01 1.35e+01 1.94e+04 -1.0 1.38e+02 - 4.77e-01 3.94e-01h 1
69 1.8680453e+01 9.45e+00 7.11e+02 -1.0 3.94e+01 - 9.63e-01 4.67e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 1.8799957e+01 1.49e+00 7.96e-02 -1.0 5.60e+01 - 1.00e+00 8.43e-01h 1
71 1.9054939e+01 2.83e-01 1.00e-03 -1.0 2.00e+01 - 1.00e+00 1.00e+00h 1
72 1.7640871e+01 1.83e+00 6.11e+05 -1.7 1.28e+01 - 9.24e-01 1.00e+00h 1
73 1.7231464e+01 3.33e-01 4.06e-03 -1.7 1.30e+01 - 1.00e+00 1.00e+00h 1
74 1.7219138e+01 1.43e-03 1.62e-06 -1.7 2.09e+00 - 1.00e+00 1.00e+00h 1
75 1.6577542e+01 6.36e-01 5.05e+05 -5.7 1.27e+01 - 7.48e-01 1.00e+00f 1
76 1.6481556e+01 8.40e-02 1.06e+05 -5.7 1.32e+01 - 7.89e-01 9.72e-01h 1
77 1.6472546e+01 1.25e-02 3.50e+03 -5.7 3.83e+00 - 9.67e-01 8.67e-01h 1
78 1.6470917e+01 1.57e-04 4.43e-08 -5.7 1.57e+00 - 1.00e+00 1.00e+00h 1
79 1.6470916e+01 1.95e-12 1.85e-11 -5.7 1.92e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.6470835e+01 3.72e-08 1.71e-10 -8.6 7.33e-03 - 1.00e+00 1.00e+00f 1
81 1.6470835e+01 9.04e-13 2.52e-14 -8.6 2.39e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 81
(scaled) (unscaled)
Objective...............: 1.6470835265947020e+01 1.6470835265947020e+01
Dual infeasibility......: 2.5209565167191177e-14 2.5209565167191177e-14
Constraint violation....: 9.0379370180735576e-13 9.0379370180735576e-13
Complementarity.........: 2.5059035896057772e-09 2.5059035896057772e-09
Overall NLP error.......: 2.5059035896057772e-09 2.5059035896057772e-09
Number of objective function evaluations = 146
Number of objective gradient evaluations = 82
Number of equality constraint evaluations = 146
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 82
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 81
Total CPU secs in IPOPT (w/o function evaluations) = 0.020
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
{'initial_tf': 0.0, 'lower_tf': 0.001, 'final_tf': 16.47083526594702, 'iterations': 79}
Solver command line: ['/Users/adowling/.idaes/bin/ipopt', '/var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmpcmb31dsh.pyomo.nl', '-AMPL']
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 769
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 135
Total number of variables............................: 228
variables with only lower bounds: 47
variables with lower and upper bounds: 91
variables with only upper bounds: 0
Total number of equality constraints.................: 184
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0999990e-02 1.00e+02 8.03e-03 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.8383675e+00 1.00e+02 1.10e+03 -1.0 9.99e+03 - 9.91e-07 1.83e-04f 1
2 3.6211597e+00 9.99e+01 5.41e+01 -1.0 3.61e+03 - 1.04e-04 4.94e-04f 1
3 4.2781033e+00 9.99e+01 5.33e+01 -1.0 5.59e+03 - 2.21e-05 1.67e-04h 1
4 4.3506853e+00 9.99e+01 5.28e+01 -1.0 2.91e+03 - 4.52e-04 3.00e-05h 1
5 4.9284606e+00 9.99e+01 5.19e+01 -1.0 3.36e+03 - 4.43e-04 2.57e-04h 1
6 5.5897468e+00 9.98e+01 5.09e+01 -1.0 3.25e+03 - 8.45e-04 3.36e-04h 2
7 6.4331918e+00 9.98e+01 4.97e+01 -1.0 3.03e+03 - 1.58e-03 4.84e-04h 1
8 6.7765220e+00 9.98e+01 5.47e+01 -1.0 2.69e+03 - 2.66e-03 2.26e-04h 3
9 6.9957257e+00 9.98e+01 6.38e+01 -1.0 2.37e+03 - 4.00e-03 1.53e-04h 4
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 7.4791752e+00 9.97e+01 7.83e+01 -1.0 2.07e+03 - 5.39e-03 3.53e-04h 4
11 7.9356105e+00 9.97e+01 1.12e+02 -1.0 1.82e+03 - 8.18e-03 3.57e-04h 5
12 8.2914405e+00 9.97e+01 1.44e+02 -1.0 1.60e+03 - 8.04e-03 2.97e-04h 6
13 8.8444541e+00 9.96e+01 1.74e+02 -1.0 1.47e+03 - 1.02e-02 4.82e-04h 6
14 9.2558537e+00 9.96e+01 1.89e+02 -1.0 1.36e+03 - 6.58e-03 3.82e-04h 7
15 9.7292439e+00 9.95e+01 2.18e+02 -1.0 1.29e+03 - 1.27e-02 4.59e-04h 7
16 1.0324823e+01 9.95e+01 2.33e+02 -1.0 1.22e+03 - 7.74e-03 6.07e-04h 7
17 1.0900785e+01 9.94e+01 2.58e+02 -1.0 1.15e+03 - 1.39e-02 6.21e-04h 7
18 5.3478817e+01 9.46e+01 7.27e+02 -1.0 1.09e+03 - 9.88e-03 4.84e-02w 1
19 5.8662871e+01 7.42e+01 2.39e+05 -1.0 1.66e+02 2.0 1.74e-04 2.15e-01w 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 9.2332994e+01 6.75e+01 1.77e+05 -1.0 3.69e+02 1.5 7.28e-05 9.13e-02w 1
21 1.1233426e+01 9.94e+01 2.75e+02 -1.0 1.07e+02 1.9 9.88e-03 3.78e-04h 7
22 1.1607544e+01 9.93e+01 3.02e+02 -1.0 1.06e+03 - 1.60e-02 4.38e-04h 8
23 1.2066316e+01 9.93e+01 3.19e+02 -1.0 1.01e+03 - 1.14e-02 5.54e-04h 8
24 1.2563656e+01 9.92e+01 3.45e+02 -1.0 9.77e+02 - 1.81e-02 6.24e-04h 8
25 1.3147773e+01 9.91e+01 3.62e+02 -1.0 9.37e+02 - 1.35e-02 7.61e-04h 8
26 1.3777638e+01 9.90e+01 3.89e+02 -1.0 8.99e+02 - 2.11e-02 8.57e-04h 8
27 1.4508302e+01 9.89e+01 4.08e+02 -1.0 9.04e+02 - 1.65e-02 1.04e-03h 8
28 1.5305844e+01 9.88e+01 4.34e+02 -1.0 9.12e+02 - 2.52e-02 1.19e-03h 8
29 1.6176968e+01 9.87e+01 4.53e+02 -1.0 9.19e+02 - 2.06e-02 1.36e-03h 8
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 1.7070337e+01 9.85e+01 4.78e+02 -1.0 9.28e+02 - 3.06e-02 1.47e-03h 8
31 9.9407952e+01 3.40e+02 1.12e+03 -1.0 9.36e+02 - 2.58e-02 1.42e-01w 1
32 9.8491768e+01 3.26e+02 1.47e+04 -1.0 1.49e+02 1.5 7.06e-04 4.02e-02w 1
33 1.0396738e+02 1.04e+01 2.02e+04 -1.0 1.34e+02 1.0 3.12e-04 9.90e-01w 1
34 9.9062113e+01 1.28e+00 1.91e+03 -1.0 6.02e+00 0.5 9.47e-01 9.90e-01h 1
35 9.8855088e+01 3.76e-03 6.27e+02 -1.0 4.48e-01 0.0 9.95e-01 1.00e+00h 1
36 9.7943077e+01 2.39e-02 3.36e-01 -1.0 9.12e-01 -0.4 1.00e+00 1.00e+00f 1
37 9.5401415e+01 1.90e-01 3.32e+05 -3.8 2.54e+00 -0.9 9.67e-01 1.00e+00f 1
38 8.8488875e+01 1.48e+00 2.81e-01 -3.8 6.91e+00 -1.4 1.00e+00 1.00e+00f 1
39 7.2625709e+01 9.06e+00 2.15e-01 -3.8 1.59e+01 -1.9 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 5.5062810e+01 1.90e+01 4.34e-01 -3.8 3.28e+01 -2.3 1.00e+00 6.94e-01f 1
41 5.2934007e+01 1.80e+01 2.62e+00 -3.8 6.13e+01 -2.8 1.00e+00 6.97e-02h 1
42 4.0901635e+01 2.41e+01 1.09e+00 -3.8 1.52e+02 -3.3 5.35e-01 2.03e-01f 1
43 4.0645374e+01 2.37e+01 4.01e+00 -3.8 4.57e+01 -2.9 1.00e+00 1.68e-02h 1
44 3.9650903e+01 1.92e+01 3.82e+00 -3.8 9.19e+01 -3.4 2.40e-02 1.83e-01h 1
45 3.0486731e+01 2.59e+01 3.86e+00 -3.8 1.22e+03 -3.8 3.96e-03 6.12e-02f 1
46 3.0448887e+01 2.58e+01 5.84e+00 -3.8 1.06e+02 -3.4 5.42e-01 3.22e-03h 1
47 3.0004789e+01 1.63e+01 5.20e+00 -3.8 4.13e+01 -3.0 1.00e+00 3.73e-01h 1
48 2.7614913e+01 1.46e+01 3.84e+00 -3.8 1.03e+02 -3.5 4.25e-02 1.82e-01h 1
49 2.5754072e+01 1.50e+01 2.73e+00 -3.8 6.11e+02 -3.9 1.26e-01 2.34e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 2.3538912e+01 1.32e+01 8.50e-01 -3.8 7.75e+01 -3.5 2.82e-01 2.06e-01h 1
51 2.2776197e+01 1.32e+01 4.76e-01 -3.8 5.54e+02 -4.0 4.61e-02 1.09e-02h 1
52 2.2490981e+01 1.28e+01 7.89e-01 -3.8 5.17e+01 -3.6 1.12e-01 3.73e-02h 1
53 2.1765815e+01 1.27e+01 7.29e-01 -3.8 4.66e+02 -4.0 2.97e-02 1.53e-02h 1
54 2.1464802e+01 1.21e+01 7.21e-01 -3.8 7.00e+01 -3.6 7.58e-02 5.27e-02h 1
55 2.0342107e+01 1.16e+01 6.87e-01 -3.8 3.80e+02 -4.1 2.83e-02 4.73e-02h 1
56 2.0131802e+01 1.07e+01 6.28e-01 -3.8 5.95e+01 -3.7 5.68e-01 7.80e-02h 1
57 1.9670710e+01 1.05e+01 6.18e-01 -3.8 3.23e+02 -4.1 9.50e-02 2.21e-02h 1
58 1.9495820e+01 9.90e+00 5.80e-01 -3.8 6.17e+01 -3.7 3.17e-01 5.96e-02h 1
59 1.8780477e+01 9.60e+00 6.03e-01 -3.8 3.25e+02 -4.2 1.22e-01 3.50e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 1.8620195e+01 8.91e+00 5.92e-01 -3.8 5.20e+01 -3.8 1.00e+00 7.25e-02h 1
61 1.8048673e+01 8.47e+00 6.01e-01 -3.8 1.62e+02 -4.2 9.69e-02 5.11e-02h 1
62 1.7885085e+01 7.46e+00 4.32e-01 -3.8 4.49e+01 -3.8 1.00e+00 1.17e-01h 1
63 1.7437290e+01 6.76e+00 4.40e-01 -3.8 1.34e+02 -4.3 1.84e-01 8.24e-02h 1
64 1.7168771e+01 6.79e+00 4.74e-01 -3.8 1.21e+03 -4.8 1.26e-02 4.59e-03h 1
65 1.6958666e+01 6.48e+00 6.41e-01 -3.8 1.17e+02 -4.3 7.77e-01 5.21e-02h 1
66 1.6634866e+01 6.46e+00 6.64e-01 -3.8 4.66e+02 -4.8 5.42e-02 1.96e-02h 1
67 1.6450158e+01 5.90e+00 5.15e-01 -3.8 7.04e+01 -4.4 1.00e+00 9.37e-02h 1
68 1.6279383e+01 5.65e+00 5.06e-01 -3.8 2.26e+02 -4.9 1.64e-01 4.84e-02h 1
69 1.6166937e+01 5.54e+00 5.26e-01 -3.8 4.19e+02 -5.3 1.35e-01 2.08e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 1.6106737e+01 5.52e+00 5.26e-01 -3.8 3.83e+03 - 5.02e-03 3.37e-03h 1
71 1.6065843e+01 5.41e+00 5.17e-01 -3.8 9.14e+02 - 4.60e-02 2.02e-02h 1
72 1.6073123e+01 5.28e+00 5.03e-01 -3.8 3.77e+02 - 4.18e-02 2.37e-02h 1
73 1.6083676e+01 4.78e+00 3.66e-01 -3.8 9.42e+01 - 1.00e+00 9.46e-02h 1
74 1.6479457e+01 1.23e-01 9.07e-03 -3.8 5.69e+00 - 1.00e+00 1.00e+00h 1
75 1.6477418e+01 1.99e-05 1.65e-06 -3.8 2.52e-01 - 1.00e+00 1.00e+00h 1
76 1.6470924e+01 2.29e-04 1.10e-06 -5.7 5.62e-01 - 1.00e+00 1.00e+00h 1
77 1.6470916e+01 6.97e-09 1.92e-11 -5.7 1.45e-02 - 1.00e+00 1.00e+00h 1
78 1.6470835e+01 3.72e-08 1.71e-10 -8.6 7.34e-03 - 1.00e+00 1.00e+00h 1
79 1.6470835e+01 1.01e-12 2.56e-14 -8.6 2.38e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 79
(scaled) (unscaled)
Objective...............: 1.6470835265947020e+01 1.6470835265947020e+01
Dual infeasibility......: 2.5565949931751082e-14 2.5565949931751082e-14
Constraint violation....: 8.9421908794747823e-13 1.0102951596350285e-12
Complementarity.........: 2.5059035896057784e-09 2.5059035896057784e-09
Overall NLP error.......: 2.5059035896057784e-09 2.5059035896057784e-09
Number of objective function evaluations = 232
Number of objective gradient evaluations = 80
Number of equality constraint evaluations = 232
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 80
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 79
Total CPU secs in IPOPT (w/o function evaluations) = 0.028
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
{'initial_tf': 1000000.0, 'lower_tf': 0.001, 'final_tf': 16.47083526718121, 'iterations': 44}
Solver command line: ['/Users/adowling/.idaes/bin/ipopt', '/var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmptleme4wu.pyomo.nl', '-AMPL']
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 769
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 135
Total number of variables............................: 228
variables with only lower bounds: 47
variables with lower and upper bounds: 91
variables with only upper bounds: 0
Total number of equality constraints.................: 184
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+06 1.00e+04 7.35e-07 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 9.9990305e+05 1.00e+04 1.17e+00 -1.0 4.48e+05 - 9.91e-05 2.17e-04f 1
2 9.9781905e+05 9.97e+03 1.02e+02 -1.0 8.12e+05 - 1.78e-04 2.57e-03f 1
3 3.7819144e+04 9.59e+03 3.61e+04 -1.0 9.98e+05 - 2.78e-03 9.62e-01f 1
4 3.6701135e+04 9.31e+03 3.50e+04 -1.0 3.78e+04 - 1.89e-02 2.96e-02f 1
5 3.6701234e+02 2.33e+03 4.70e+02 -1.0 3.67e+04 - 5.09e-02 9.90e-01f 1
6 7.5310045e+01 1.97e+03 2.10e+04 -1.0 5.69e+02 - 1.20e-02 8.08e-01f 1
7 8.4388532e+01 1.85e+03 1.66e+04 -1.0 8.16e+02 - 8.17e-02 5.79e-02h 1
8 9.0827033e+01 1.78e+03 1.55e+04 -1.0 8.78e+02 - 3.78e-02 3.41e-02h 1
9 1.0713654e+02 1.53e+03 2.85e+04 -1.0 7.51e+02 - 3.16e-02 1.33e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.1481720e+02 1.44e+03 1.20e+04 -1.0 4.17e+02 - 1.45e-01 5.87e-02h 1
11 1.3668225e+02 1.11e+03 2.76e+04 -1.0 5.38e+02 - 1.21e-01 2.21e-01h 1
12 1.2429708e+02 5.71e+02 1.13e+05 -1.0 1.31e+03 - 5.67e-02 4.86e-01f 1
13 1.2442612e+02 5.49e+02 3.56e+03 -1.0 5.27e+02 - 2.78e-01 3.83e-02h 1
14 1.2745970e+02 4.52e+02 7.21e+03 -1.0 4.01e+02 - 2.00e-01 1.77e-01h 1
15 1.0761486e+02 3.40e+02 1.18e+05 -1.0 1.27e+03 - 3.22e-03 2.34e-01f 1
16 1.0933068e+02 3.04e+02 2.84e+04 -1.0 2.12e+02 - 2.19e-01 1.07e-01h 1
17 9.9340849e+01 2.53e+02 1.43e+05 -1.0 9.85e+02 - 4.66e-03 1.62e-01f 1
18 1.1172284e+02 5.74e+00 6.64e+05 -1.0 1.60e+02 - 2.50e-01 1.00e+00h 1
19 8.6975271e+01 6.13e+01 3.74e+06 -1.0 4.38e+02 - 1.80e-01 5.49e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 8.6217009e+01 5.74e+01 2.72e+06 -1.0 3.33e+02 - 1.66e-01 6.40e-02f 1
21 8.9317885e+01 1.17e+01 2.84e+06 -1.0 7.10e+01 - 4.46e-01 8.11e-01f 1
22 8.5320799e+01 1.14e+01 2.66e+06 -1.0 2.98e+02 - 9.36e-02 2.14e-01f 1
23 7.9309753e+01 1.42e+01 1.95e+06 -1.0 5.07e+02 - 2.50e-01 1.82e-01f 1
24 7.3201186e+01 1.66e+01 1.71e+06 -1.0 5.99e+02 - 1.33e-01 1.70e-01f 1
25 6.5845444e+01 1.58e+01 1.65e+06 -1.0 2.64e+02 - 1.02e-01 4.13e-01f 1
26 6.1736897e+01 1.37e+01 1.50e+06 -1.0 3.06e+02 - 1.13e-01 2.61e-01f 1
27 5.8715499e+01 1.36e+01 1.22e+06 -1.0 3.60e+02 - 1.76e-01 1.13e-01f 1
28 5.4245993e+01 1.50e+01 1.11e+06 -1.0 5.06e+02 - 9.79e-02 1.45e-01f 1
29 4.9891552e+01 1.36e+01 9.68e+05 -1.0 2.14e+02 - 1.43e-01 2.87e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 4.5165997e+01 1.21e+01 8.11e+05 -1.0 2.78e+02 - 1.82e-01 3.96e-01f 1
31 3.8709430e+01 1.38e+01 7.04e+05 -1.0 2.22e+02 - 1.56e-01 5.34e-01f 1
32 3.5158500e+01 1.32e+01 5.31e+05 -1.0 1.84e+02 - 2.46e-01 2.92e-01f 1
33 3.0165631e+01 1.45e+01 3.85e+05 -1.0 1.76e+02 - 2.79e-01 3.96e-01f 1
34 2.4558035e+01 1.62e+01 2.57e+05 -1.0 1.32e+02 - 3.36e-01 5.24e-01f 1
35 1.9946467e+01 1.67e+01 1.57e+05 -1.0 1.02e+02 - 3.92e-01 5.59e-01f 1
36 1.8355638e+01 1.31e+01 3.16e+04 -1.0 7.07e+01 - 7.93e-01 3.37e-01f 1
37 1.8671537e+01 6.08e-01 4.08e-02 -1.0 3.36e+01 - 1.00e+00 1.00e+00h 1
38 1.7216736e+01 2.32e+00 3.24e+06 -5.7 1.93e+01 - 6.76e-01 1.00e+00f 1
39 1.6791543e+01 8.58e-01 1.09e+06 -5.7 1.87e+01 - 6.63e-01 7.89e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 1.6580166e+01 3.21e-01 3.04e+05 -5.7 1.74e+01 - 7.21e-01 7.65e-01f 1
41 1.6483937e+01 1.30e-01 4.91e+04 -5.7 2.17e+01 - 8.39e-01 9.24e-01f 1
42 1.6471275e+01 3.97e-03 1.57e+03 -5.7 2.32e+00 - 9.68e-01 9.69e-01f 1
43 1.6470916e+01 4.84e-07 8.17e-09 -5.7 1.82e-02 - 1.00e+00 1.00e+00h 1
44 1.6470835e+01 3.70e-08 1.70e-10 -8.6 7.32e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 44
(scaled) (unscaled)
Objective...............: 1.6470835267181211e+01 1.6470835267181211e+01
Dual infeasibility......: 1.6982534595702136e-10 1.6982534595702136e-10
Constraint violation....: 3.6981354867293698e-12 3.6981354867293703e-08
Complementarity.........: 3.2032414677344006e-09 3.2032414677344006e-09
Overall NLP error.......: 3.2032414677344006e-09 3.6981354867293703e-08
Number of objective function evaluations = 45
Number of objective gradient evaluations = 45
Number of equality constraint evaluations = 45
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 45
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 44
Total CPU secs in IPOPT (w/o function evaluations) = 0.011
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
Another Version¶
Another modeling option is to define an extra differential equation to compute time. The normalized-time formulation is unchanged: tau remains the fixed independent variable on , and tf remains a decision variable. The additional state is useful when physical time must be available along the trajectory; otherwise, recover it after solving as .
def create_model2():
# Define the model
m = pyo.ConcreteModel()
# Define the model parameters
m.R = pyo.Param(initialize=0.001, units=1 / units.m) # Friction factor
m.L = pyo.Param(initialize=100.0, units=units.m) # Final position
# Define time
m.tau = dae.ContinuousSet(bounds=(0, 1)) # Dimensionless time set
m.time = pyo.Var(m.tau, bounds=(0, None), units=units.s) # Time
m.tf = pyo.Var(initialize=1, bounds=(1e-3, None), units=units.s) # Final time
# Define remaining algebraic variables
m.x = pyo.Var(m.tau, bounds=(0, m.L + 50 * units.m), units=units.m) # Position
m.v = pyo.Var(m.tau, bounds=(0, None), units=units.m / units.s) # Velocity
m.u = pyo.Var(
m.tau, bounds=(-3.0, 1.0), initialize=0, units=units.m / units.s / units.s
) # Acceleration
# Define derivative variables. tau is dimensionless, so d/dtau carries the
# same units as the state. DerivativeVar does not infer these -- say them.
m.dtime = dae.DerivativeVar(m.time, units=units.s)
m.dx = dae.DerivativeVar(m.x, units=units.m)
m.dv = dae.DerivativeVar(m.v, units=units.m / units.s)
# Declare the objective (minimize final time)
m.obj = pyo.Objective(expr=m.tf)
# Define the constraints
# position
@m.Constraint(m.tau)
def ode1(m, i):
if i == 0:
return pyo.Constraint.Skip
return m.dx[i] == m.tf * m.v[i]
# velocity
@m.Constraint(m.tau)
def ode2(m, i):
if i == 0:
return pyo.Constraint.Skip
return m.dv[i] == m.tf * (m.u[i] - m.R * m.v[i] ** 2)
# time
@m.Constraint(m.tau)
def ode3(m, i):
if i == 0:
return pyo.Constraint.Skip
return m.dtime[i] == m.tf
# Define the initial/boundary conditions
def _init(m):
yield m.x[0] == 0
yield m.x[1] == m.L
yield m.v[0] == 0
yield m.v[1] == 0
yield m.time[0] == 0
m.initcon = pyo.ConstraintList(rule=_init)
# Units are checked on the continuous model, before discretization.
assert_units_consistent(m)
return m
m = create_model2()# Declare the discretizer
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m, nfe=15, scheme="LAGRANGE-RADAU", ncp=3)
# force piecewise constant controls (acceleration) over each finite element
m = discretizer.reduce_collocation_points(m, var=m.u, ncp=1, contset=m.tau)
# Solve
solver = pyo.SolverFactory("ipopt")
results = solver.solve(m, tee=True)
assert pyo.check_optimal_termination(results), (
f"Solve failed: status={results.solver.status}, "
f"termination={results.solver.termination_condition}"
)
print("final time = %6.2f seconds" % (pyo.value(m.tf)))Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1145
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 135
Total number of variables............................: 319
variables with only lower bounds: 93
variables with lower and upper bounds: 91
variables with only upper bounds: 0
Total number of equality constraints.................: 305
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 1.00e+02 3.64e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 9.9963150e-01 9.99e+01 1.09e+02 -1.0 1.00e+03 - 7.69e-05 1.26e-03f 1
2 1.0670594e+00 9.98e+01 1.21e+02 -1.0 8.13e+02 - 1.78e-04 1.13e-03f 1
3 1.6605461e+00 9.94e+01 3.87e+01 -1.0 6.32e+02 - 4.66e-04 3.47e-03f 1
4 2.1544627e+00 9.89e+01 2.22e+01 -1.0 3.71e+02 - 5.93e-03 4.72e-03f 1
5 3.3655947e+00 9.75e+01 2.07e+01 -1.0 2.82e+02 - 6.66e-03 1.47e-02f 1
6 3.5350678e+00 9.74e+01 2.03e+01 -1.0 3.04e+03 -2.0 1.12e-04 5.44e-04f 3
7 6.5718493e+00 9.24e+01 4.59e+01 -1.0 2.00e+02 -0.7 3.31e-02 5.18e-02f 1
8 2.0354039e+01 8.38e+01 3.17e+02 -1.0 2.04e+02 -1.1 8.88e-02 4.11e-01f 1
9 2.0635680e+01 7.80e+01 2.96e+02 -1.0 1.35e+02 -1.6 3.64e-02 6.86e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.7964072e+01 3.94e+01 2.80e+02 -1.0 3.06e+02 -1.2 3.55e-03 2.25e-01f 1
11 3.3187087e+01 8.22e+00 8.36e+01 -1.0 9.77e+01 -1.7 5.82e-02 9.90e-01h 1
12 3.3005960e+01 1.18e-01 1.27e+00 -1.0 2.61e+00 -2.2 9.91e-01 9.90e-01h 1
13 2.1889688e+01 3.95e+01 2.64e+04 -1.0 4.08e+02 - 2.58e-01 3.61e-01f 1
14 2.0331592e+01 2.48e+01 1.17e+05 -1.0 6.35e+01 - 7.55e-01 4.17e-01h 1
15 1.9469251e+01 3.80e+00 1.13e+05 -1.0 9.54e+01 - 8.22e-01 1.00e+00h 1
16 2.0063941e+01 1.36e+00 1.84e-02 -1.0 4.30e+01 - 1.00e+00 1.00e+00h 1
17 1.9748189e+01 2.25e-01 9.31e-04 -1.0 1.16e+01 - 1.00e+00 1.00e+00h 1
18 1.7573422e+01 3.82e+00 8.47e+05 -1.7 1.75e+01 - 8.94e-01 1.00e+00h 1
19 1.7308137e+01 1.91e-01 3.38e-03 -1.7 8.97e+00 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.6682427e+01 5.30e-01 8.81e+04 -2.5 7.99e+00 - 9.49e-01 1.00e+00h 1
21 1.6636675e+01 1.02e-02 4.20e-05 -2.5 3.51e+00 - 1.00e+00 1.00e+00h 1
22 1.6527711e+01 1.85e-02 3.25e-04 -3.8 1.58e+00 - 1.00e+00 1.00e+00h 1
23 1.6526505e+01 2.59e-06 3.09e-08 -3.8 2.42e-02 - 1.00e+00 1.00e+00h 1
24 1.6520273e+01 6.33e-05 1.05e-06 -5.7 9.59e-02 - 1.00e+00 1.00e+00h 1
25 1.6520269e+01 3.73e-11 1.87e-11 -5.7 1.11e-04 - 1.00e+00 1.00e+00h 1
26 1.6520192e+01 9.76e-09 1.61e-10 -8.6 1.19e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 26
(scaled) (unscaled)
Objective...............: 1.6520191520684008e+01 1.6520191520684008e+01
Dual infeasibility......: 1.6135454785868610e-10 1.6135454785868610e-10
Constraint violation....: 9.7605266091704834e-09 9.7605266091704834e-09
Complementarity.........: 2.5758043844751979e-09 2.5758043844751979e-09
Overall NLP error.......: 9.7605266091704834e-09 9.7605266091704834e-09
Number of objective function evaluations = 30
Number of objective gradient evaluations = 27
Number of equality constraint evaluations = 30
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 27
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 26
Total CPU secs in IPOPT (w/o function evaluations) = 0.012
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
final time = 16.52 seconds