Dynamic Optimization with Pyomo.DAE: An Example
Contents
3.4. Dynamic Optimization with Pyomo.DAE: An Example#
Prepared by: Prof. Alexander Dowling, Molly Dougher (mdoughe6@nd.edu, 2023)
import sys
if "google.colab" in sys.modules:
!wget "https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py"
import helper
#!pip install casadi
helper.install_idaes()
helper.install_ipopt()
--2023-05-04 21:36:24-- https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 7171 (7.0K) [text/plain]
Saving to: ‘helper.py’
helper.py 0%[ ] 0 --.-KB/s
helper.py 100%[===================>] 7.00K --.-KB/s in 0s
2023-05-04 21:36:24 (53.9 MB/s) - ‘helper.py’ saved [7171/7171]
Installing idaes via pip...
idaes was successfully installed
Running idaes get-extensions to install Ipopt, k_aug, and more...
Ipopt 3.13.2 (x86_64-pc-linux-gnu), ASL(20190605)
[K_AUG] 0.1.0, Part of the IDAES PSE framework
Please visit https://idaes.org/ (x86_64-pc-linux-gnu), ASL(20190605)
Couenne 0.5.8 -- an Open-Source solver for Mixed Integer Nonlinear Optimization
Mailing list: couenne@list.coin-or.org
Instructions: http://www.coin-or.org/Couenne
couenne (x86_64-pc-linux-gnu), ASL(20190605)
Bonmin 1.8.8 using Cbc 2.10.8 and Ipopt 3.13.2
bonmin (x86_64-pc-linux-gnu), ASL(20190605)
Ipopt 3.13.3 (x86_64-pc-linux-gnu), ASL(20190605)
ipopt was successfully installed
k_aug was successfuly installed
cbc was successfuly installed
clp was successfuly installed
bonmin was successfuly installed
couenne was successfuly installed
ipopt_l1 was successfuly installed
3.4.1. Car Example#
Adapted from https://github.com/Pyomo/pyomo/blob/master/examples/dae/car_example.py
You are a race car driver with a simple goal. Drive distance \(L\) in the minimal amount of time but come to a complete stop at the finish line.
3.4.1.1. Optimal Control Problem Formulation#
Mathematically, you want to solve the following optimal control problem:
\[\begin{split}\begin{align*}
\min_{u} \quad & t_f \\
\mathrm{s.t.} \quad & \frac{dx}{dt} = v \\
& \frac{dv}{dt} = u - R v^2 \\
& x(t=0) = 0, ~~ x(t=t_f) = L \\
& v(t=0) = 0, ~~ v(t=t_f) = 0 \\
& -3 \leq u \leq 1
\end{align*}\end{split}\]
where \(a\) is the acceleration/braking (your control variable) and \(R\) is the drag coefficient (parameter).
3.4.1.2. Scale Time#
Let \(t = \tau \cdot t_f\) where \(\tau \in [0,1]\). Thus \(dt = t_f d\tau\). The optimal control problem becomes:
\[\begin{split}\begin{align*}
\min_{u} \quad & t_f \\
\mathrm{s.t.} \quad & \frac{dx}{d\tau} = t_f v \\
& \frac{dv}{d\tau} = t_f (u - R v^2) \\
& x(\tau = 0) = 0, ~~ x(\tau = 1) = L \\
& v(\tau = 0) = 0, ~~ v(\tau = 1) = 0 \\
& -3 \leq u \leq 1
\end{align*}\end{split}\]
3.4.1.3. Orthogonal Collocation on Finite Elements: Manual Approach#
Here is the “classical”/”old school” way of manually implementing collocation on finite elements with sets.
Indices and Sets
Finite elements: \(i \in \mathcal{I} = \{1,2,...,N_{FE}\}\)
Finite elements (except first): \(\mathcal{I} \setminus 1\)
Internal collocation points: \(j,k \in \mathcal{J} = \{1,2,...,N_{C}\}\)
Parameters
Coefficients in collocation/Runge-Kutta formula: \(a_{j,j}\)
Drag coefficient: \(R\)
Race length: \(L\)
Scaled time for each finite element: \(h_i = \frac{1}{N_{FE}}\)
Variables
Position (internal collocation points): \(x_{i,j}\)
Position (beginning of each finite element): \(\bar{x}_{i}\)
Position derivative (internal collocation points): \(\dot{x}_{i,j}\)
Velocity (internal collocation points): \(v_{i,j}\)
Velocity (beginning of each finite element): \(\bar{v}_{i}\)
Velocity derivative (internal collocation points): \(\dot{v}_{i,j}\)
Time (internal collocation points): \(t_{i,j}\)
Time (beginning of each finite element): \(\bar{t}_{i}\)
Acceleration (control degrees of freedom): \(u_i\)
Objective and Constraints
\[\begin{split}\begin{align*}
\min \quad & t_f \\
\mathrm{s.t.} \quad & \mathrm{Differential~equation~for~x:} \\
& x_{i,j} = \bar{x}_{i} + h_i \sum_{k \in \mathcal{J}} a_{k,j} \dot{x}_{i,k},\quad \forall i \in \mathcal{I},~j \in \mathcal{J} \\
& \bar{x}_{i} = \bar{x}_{i-1} + h_{i-1} \sum_{k \in \mathcal{J}} a_{k,N_C} \dot{x}_{i-1,k},\quad \forall i \in \mathcal{I} \setminus 1\\
& \dot{x}_{i,j} = t_f v_{i,j}, \quad \forall i \in \mathcal{I} \setminus 1,~j \in \mathcal{J} \\ \\
& \mathrm{Differential~equation~for~v:} \\
& v_{i,j} = \bar{v}_{i} + h_i \sum_{k \in \mathcal{J}} a_{k,j} \dot{v}_{i,k},\quad \forall i \in \mathcal{I},~j \in \mathcal{J} \\
& \bar{v}_{i} = \bar{v}_{i-1} + h_{i-1} \sum_{k \in \mathcal{J}} a_{k,N_C} \dot{v}_{i-1,k},\quad \forall i \in \mathcal{I} \setminus 1\\
& \dot{v}_{i,j} = t_f (u_{i} - R v_{i,j}^2), \quad \forall i \in \mathcal{I} \setminus 1,~j \in \mathcal{J} \\ \\
& \mathrm{Differential~equation~for~time:} \\
& t_{i,j} = \bar{t}_{i} + h_i \sum_{k \in \mathcal{J}} a_{k,j} t_f,\quad \forall i \in \mathcal{I},~j \in \mathcal{J} \\
& \bar{t}_{i} = \bar{t}_{i-1} + h_{i-1} \sum_{k \in \mathcal{J}} a_{k,N_C} t_f,\quad \forall i \in \mathcal{I} \setminus 1\\ \\
& \mathrm{Initial~conditions:} \\
& \bar{x}_{1} = 0, \quad \bar{v}_{1} = 0, \quad \bar{t}_1 = 0\\ \\
& \mathrm{Boundary~conditions:} \\
& L = \bar{x}_{N_{FE}} + h_{N_{FE}} \sum_{k \in \mathcal{J}} a_{k,N_C} \dot{x}_{N_{FE},k}, \\
& 0 = \bar{v}_{N_{FE}} + h_{N_{FE}} \sum_{k \in \mathcal{J}} a_{k,N_C} \dot{v}_{N_{FE},k}, \\ \\
& \mathrm{Bounds~on~acceleration:} \\
& -3 \leq a_i \leq 1, \quad \forall i \in \mathcal{I}
\end{align*}\end{split}\]
import pyomo.environ as pyo
import matplotlib.pyplot as plt
import pyomo.dae as dae
# Define the model
m2 = pyo.ConcreteModel()
# Define model parameters
m2.R = pyo.Param(initialize=0.001) # Friction factor
m2.L = pyo.Param(initialize=100.0) # Final position
# Define finite elements and collocation points
NFE = 15 # Number of finite elements
NC = 3 # Number of collocation points
m2.I = pyo.Set(initialize=pyo.RangeSet(1,NFE)) # Set of finite elements
m2.J = pyo.Set(initialize=pyo.RangeSet(1,NC)) # Set of internal collocation points
# Define first order derivative collocation matrix
A = {}
A[1,1] = 0.19681547722366
A[1,2] = 0.39442431473909
A[1,3] = 0.37640306270047
A[2,1] = -0.06553542585020
A[2,2] = 0.29207341166523
A[2,3] = 0.51248582618842
A[3,1] = 0.02377097434822
A[3,2] = -0.04154875212600
A[3,3] = 0.11111111111111
# Define A matrix as a model parameter
m2.a = pyo.Param(m2.J, m2.J, initialize=A)
# Define step for each finite element
m2.h = pyo.Param(m2.I,initialize=1/NFE)
# Define objective (final time)
m2.tf = pyo.Var(domain=pyo.NonNegativeReals)
# Variables for x (position)
m2.x0 = pyo.Var(m2.I) # Internal collocation points
m2.x = pyo.Var(m2.I,m2.J) # Beginning of each finite element
m2.xdot = pyo.Var(m2.I,m2.J) # Derivative, beginning of each finite element
# Variables for v (velocity)
m2.v0 = pyo.Var(m2.I) # Internal collocation points
m2.v = pyo.Var(m2.I,m2.J) # Beginning of each finite element
m2.vdot = pyo.Var(m2.I,m2.J) # Derivative, beginning of each finite element
# Variables for t
m2.t0 = pyo.Var(m2.I) # Internal collocation points
m2.t = pyo.Var(m2.I,m2.J) # Beginning of each finite element
# Acceleration
m2.u = pyo.Var(m2.I, bounds=(-3,1)) # Control DOF
### Finite Element Collocation Equations
# position
def FECOLx_(m2,i,j):
return m2.x[i,j] == m2.x0[i] + m2.h[i]*sum(m2.a[k,j]*m2.xdot[i,k] for k in m2.J)
m2.FECOLx = pyo.Constraint(m2.I,m2.J,rule=FECOLx_)
# velocity
def FECOLv_(m2,i,j):
return m2.v[i,j] == m2.v0[i] + m2.h[i]*sum(m2.a[k,j]*m2.vdot[i,k] for k in m2.J)
m2.FECOLv = pyo.Constraint(m2.I,m2.J,rule=FECOLv_)
# time
def FECOLt_(m2,i,j):
return m2.t[i,j] == m2.t0[i] + m2.h[i]*sum(m2.a[k,j]*m2.tf for k in m2.J)
m2.FECOLt = pyo.Constraint(m2.I,m2.J,rule=FECOLt_)
### Continuity Equations
# position
def CONx_(m2,i):
if i == 1:
return pyo.Constraint.Skip
else:
return m2.x0[i] == m2.x0[i-1] + m2.h[i-1]*sum(m2.a[k,NC]*m2.xdot[i,k] for k in m2.J)
m2.CONx = pyo.Constraint(m2.I, rule=CONx_)
# velocity
def CONv_(m2,i):
if i == 1:
return pyo.Constraint.Skip
else:
return m2.v0[i] == m2.v0[i-1] + m2.h[i-1]*sum(m2.a[k,NC]*m2.vdot[i,k] for k in m2.J)
m2.CONv = pyo.Constraint(m2.I, rule=CONv_)
# time
def CONt_(m2,i):
if i == 1:
return pyo.Constraint.Skip
else:
return m2.t0[i] == m2.t0[i-1] + m2.h[i-1]*sum(m2.a[k,NC]*m2.tf for k in m2.J)
m2.CONt = pyo.Constraint(m2.I, rule=CONt_)
### Differential equations
# position
def ODEx_(m2,i,j):
return m2.xdot[i,j] == m2.tf*m2.v[i,j]
m2.ODEx = pyo.Constraint(m2.I,m2.J,rule=ODEx_)
# velocity
def ODEv_(m2,i,j):
return m2.vdot[i,j] == m2.tf*(m2.u[i] - m2.R*m2.v[i,j]**2)
m2.ODEv = pyo.Constraint(m2.I,m2.J,rule=ODEv_)
### Initial conditions
m2.xIC = pyo.Constraint(expr=m2.x0[1] == 0)
m2.vIC = pyo.Constraint(expr=m2.v0[1] == 0)
m2.tIC = pyo.Constraint(expr=m2.t0[1] == 0)
### Boundary conditions
m2.xBC = pyo.Constraint(expr=m2.L == m2.x0[NFE] + m2.h[NFE]*sum(m2.a[k,NC]*m2.xdot[NFE,k] for k in m2.J))
m2.vBC = pyo.Constraint(expr=0 == m2.v0[NFE] + m2.h[NFE]*sum(m2.a[k,NC]*m2.vdot[NFE,k] for k in m2.J))
### Set objective
m2.obj = pyo.Objective(expr=m2.tf)
# Solve the model
solver = pyo.SolverFactory('ipopt')
solver.solve(m2,tee=True)
print("final time = %6.2f" %(pyo.value(m2.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...: 1093
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 105
Total number of variables............................: 286
variables with only lower bounds: 1
variables with lower and upper bounds: 15
variables with only upper bounds: 0
Total number of equality constraints.................: 272
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 9.9999900e-03 1.00e+02 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1r 9.9999900e-03 1.00e+02 9.99e+02 2.0 0.00e+00 - 0.00e+00 1.05e-07R 2
2r 3.0578540e+00 9.99e+01 9.31e+02 2.0 4.49e+01 - 1.35e-01 6.80e-02f 2
3r 2.2146807e+00 9.91e+01 5.68e+02 1.3 2.16e+00 - 1.00e+00 3.91e-01f 1
4r 2.0626784e+00 9.87e+01 4.21e+02 0.6 2.94e+00 - 6.36e-01 2.59e-01f 1
5r 2.7501672e+00 9.77e+01 2.42e+02 0.6 8.74e+00 - 5.56e-01 4.24e-01f 1
6r 5.8696920e+00 9.43e+01 1.91e+02 -0.1 8.62e+00 - 4.13e-01 7.57e-01f 1
7r 8.6536528e+00 8.80e+01 1.65e+02 -0.1 1.29e+01 - 2.98e-01 7.27e-01f 1
8r 9.0137628e+00 8.55e+01 1.09e+02 -0.1 3.77e+00 - 5.42e-01 1.00e+00f 1
9 9.4134225e+00 8.42e+01 9.84e-01 -1.0 1.13e+02 - 3.63e-01 1.56e-02h 7
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 9.7719760e+00 8.29e+01 9.61e-01 -1.0 1.13e+02 - 6.94e-01 1.56e-02h 7
11 1.0069154e+01 8.16e+01 1.10e+00 -1.0 1.20e+02 - 6.26e-01 1.52e-02h 7
12 1.0350049e+01 8.04e+01 1.27e+00 -1.0 1.24e+02 - 1.00e+00 1.56e-02h 7
13 1.0613260e+01 7.91e+01 1.29e+00 -1.0 1.25e+02 - 6.03e-01 1.56e-02h 7
14 1.0858870e+01 7.79e+01 1.31e+00 -1.0 1.25e+02 - 1.00e+00 1.56e-02h 7
15 1.1089714e+01 7.66e+01 1.31e+00 -1.0 1.23e+02 - 7.24e-01 1.56e-02h 7
16 1.1306240e+01 7.55e+01 1.31e+00 -1.0 1.24e+02 - 1.00e+00 1.56e-02h 7
17 1.1714847e+01 7.31e+01 1.28e+00 -1.0 1.22e+02 - 8.51e-01 3.12e-02h 6
18 1.2077191e+01 7.08e+01 1.24e+00 -1.0 1.23e+02 - 1.00e+00 3.12e-02h 6
19 2.2531170e+01 7.10e+01 8.96e-01 -1.0 1.19e+02 - 1.00e+00 1.00e+00w 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.7706786e+01 4.49e+00 1.75e-01 -1.0 2.28e+01 - 1.00e+00 1.00e+00w 1
21 1.7349863e+01 1.01e-01 7.40e-04 -1.0 8.28e+00 - 1.00e+00 1.00e+00h 1
22 1.6478946e+01 1.06e+00 3.16e-02 -2.5 1.68e+01 - 8.21e-01 1.00e+00h 1
23 1.6333670e+01 1.11e-01 9.02e-04 -2.5 1.27e+01 - 9.80e-01 1.00e+00h 1
24 1.6324255e+01 5.10e-04 4.40e-06 -2.5 6.64e-01 - 1.00e+00 1.00e+00h 1
25 1.6289923e+01 3.05e-03 1.05e-05 -3.8 1.03e+00 - 1.00e+00 1.00e+00h 1
26 1.6289753e+01 1.16e-07 1.47e-09 -3.8 5.97e-03 - 1.00e+00 1.00e+00h 1
27 1.6287822e+01 9.84e-06 3.34e-08 -5.7 5.95e-02 - 1.00e+00 1.00e+00f 1
28 1.6287798e+01 1.59e-09 5.42e-12 -8.6 7.53e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 28
(scaled) (unscaled)
Objective...............: 1.6287797665173553e+01 1.6287797665173553e+01
Dual infeasibility......: 5.4232879425554638e-12 5.4232879425554638e-12
Constraint violation....: 1.5949126463965513e-09 1.5949126463965513e-09
Complementarity.........: 2.5440456242856559e-09 2.5440456242856559e-09
Overall NLP error.......: 2.5440456242856559e-09 2.5440456242856559e-09
Number of objective function evaluations = 106
Number of objective gradient evaluations = 23
Number of equality constraint evaluations = 106
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 30
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 28
Total CPU secs in IPOPT (w/o function evaluations) = 0.040
Total CPU secs in NLP function evaluations = 0.003
EXIT: Optimal Solution Found.
final time = 16.29
3.4.1.4. Orthogonal Collocation on Finite Elements: Pyomo.dae#
There is a better way! We can use Pyomo.dae to automatically formulate the collocation equations.
3.4.1.4.1. Declare Model#
# Define the model
m = pyo.ConcreteModel()
# Deine the model parameters
m.R = pyo.Param(initialize=0.001) # Friction factor
m.L = pyo.Param(initialize=100.0) # Final position
# Define time
m.tau = dae.ContinuousSet(bounds=(0,1)) # Unscaled time
m.time = pyo.Var(m.tau) # Scaled time
m.tf = pyo.Var() # Final time
# Define remaining algebraic variables
m.x = pyo.Var(m.tau,bounds=(0,m.L+50)) # Position
m.v = pyo.Var(m.tau,bounds=(0,None)) # Velocity
m.u = pyo.Var(m.tau, bounds=(-3.0,1.0),initialize=0) # Acceleration
# Define derivative variables
m.dtime = dae.DerivativeVar(m.time)
m.dx = dae.DerivativeVar(m.x)
m.dv = dae.DerivativeVar(m.v)
# Declare the objective (minimize final time)
m.obj = pyo.Objective(expr=m.tf)
# Define the constraints
# position
def _ode1(m,i):
if i == 0 :
return pyo.Constraint.Skip
return m.dx[i] == m.tf * m.v[i]
m.ode1 = pyo.Constraint(m.tau, rule=_ode1)
# velocity
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)
m.ode2 = pyo.Constraint(m.tau, rule=_ode2)
# time
def _ode3(m,i):
if i == 0:
return pyo.Constraint.Skip
return m.dtime[i] == m.tf
m.ode3 = pyo.Constraint(m.tau, rule=_ode3)
# Define the inital/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)
3.4.1.4.2. Discretize/Transcribe and Solve#
#discretizer = TransformationFactory('dae.finite_difference')
#discretizer.apply_to(m,nfe=15,scheme='BACKWARD')
# Declare the discretizer
discretizer = pyo.TransformationFactory('dae.collocation')
discretizer.apply_to(m,nfe=15,scheme='LAGRANGE-RADAU',ncp=3)
#discretizer.apply_to(m,nfe=15,scheme='LAGRANGE-LEGENDRE',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')
solver.solve(m,tee=True)
print("final time = %6.2f" %(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: 46
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 0.0000000e+00 1.00e+02 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.4627512e+01 9.96e+01 6.56e+12 -1.0 1.00e+04 - 9.91e-05 3.46e-03h 9
2 3.4471202e+01 8.66e+00 5.97e+11 -1.0 1.29e+02 8.0 2.71e-03 9.13e-01f 1
3 3.4467989e+01 7.88e+00 5.43e+11 -1.0 1.11e+01 8.4 9.15e-01 9.05e-02h 1
4 3.4467121e+01 7.88e+00 5.43e+11 -1.0 1.02e+01 7.9 9.98e-01 5.69e-04h 1
5 3.4467121e+01 7.88e+00 5.46e+11 -1.0 1.01e+01 8.4 1.00e+00 1.03e-05h 1
6 3.4504247e+01 7.30e+00 5.44e+11 -1.0 1.08e+01 7.9 1.00e+00 7.37e-02h 1
7 3.5054685e+01 2.76e+00 2.10e+11 -1.0 9.99e+00 8.3 1.00e+00 6.22e-01h 1
8 3.5393144e+01 2.75e-02 2.11e+09 -1.0 3.84e+00 7.8 1.00e+00 9.90e-01h 1
9 3.5395762e+01 1.76e-04 2.14e+07 -1.0 4.12e-02 7.4 9.92e-01 9.94e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 3.5395779e+01 2.26e-11 9.25e+05 -1.0 2.63e-04 6.9 9.90e-01 1.00e+00f 1
11 3.5395780e+01 1.82e-12 5.65e+04 -1.0 4.00e-06 6.4 9.91e-01 1.00e+00f 1
12 3.5395875e+01 1.62e-09 3.33e+02 -1.0 3.84e-04 5.9 1.00e+00 1.00e+00f 1
13 3.5395879e+01 2.80e-12 4.64e+00 -1.7 1.61e-05 5.5 1.00e+00 1.00e+00f 1
14 3.5395879e+01 1.82e-12 1.25e-01 -3.8 1.30e-06 5.0 1.00e+00 1.00e+00h 1
15 3.5395879e+01 2.73e-12 1.35e-01 -3.8 4.20e-06 4.5 1.00e+00 1.00e+00h 1
16 3.5395881e+01 2.73e-12 1.32e-01 -3.8 1.23e-05 4.0 1.00e+00 1.00e+00h 1
17 3.5395885e+01 5.49e-12 1.24e-01 -3.8 3.48e-05 3.6 1.00e+00 1.00e+00f 1
18 3.5395895e+01 3.04e-11 1.07e-01 -3.8 9.04e-05 3.1 1.00e+00 1.00e+00f 1
19 3.5395909e+01 9.38e-11 8.34e-02 -3.8 2.10e-04 2.6 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 3.5395908e+01 5.45e-12 6.13e-02 -3.8 4.64e-04 2.1 1.00e+00 1.00e+00f 1
21 3.5395824e+01 2.03e-09 4.63e-02 -3.8 1.05e-03 1.6 1.00e+00 1.00e+00f 1
22 3.5395429e+01 1.84e-08 3.79e-02 -3.8 2.58e-03 1.2 1.00e+00 1.00e+00f 1
23 3.5394058e+01 2.13e-07 3.41e-02 -3.8 6.97e-03 0.7 1.00e+00 1.00e+00f 1
24 3.5389786e+01 2.06e-06 3.31e-02 -3.8 2.03e-02 0.2 1.00e+00 1.00e+00f 1
25 3.5376906e+01 1.87e-05 3.29e-02 -3.8 6.06e-02 -0.3 1.00e+00 1.00e+00f 1
26 3.5338245e+01 1.69e-04 3.29e-02 -3.8 1.82e-01 -0.7 1.00e+00 1.00e+00f 1
27 3.5222163e+01 1.52e-03 3.29e-02 -3.8 5.45e-01 -1.2 1.00e+00 1.00e+00f 1
28 3.4873033e+01 1.39e-02 3.29e-02 -3.8 1.63e+00 -1.7 1.00e+00 1.00e+00f 1
29 3.3817946e+01 1.29e-01 3.28e-02 -3.8 4.88e+00 -2.2 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 3.0591122e+01 1.29e+00 3.23e-02 -3.8 1.44e+01 -2.7 1.00e+00 1.00e+00f 1
31 2.5066158e+01 5.19e+00 7.20e-01 -3.8 3.97e+01 -3.1 1.00e+00 5.57e-01f 1
32 2.3501854e+01 5.15e+00 1.16e+00 -3.8 9.48e+01 -3.6 1.00e+00 1.28e-01h 1
33 2.1055819e+01 6.43e+00 1.29e+00 -3.8 2.18e+02 -4.1 2.34e-01 1.16e-01h 1
34 1.9862481e+01 4.65e+00 3.98e-01 -3.8 4.62e+01 -3.7 1.00e+00 4.13e-01h 1
35 1.8481817e+01 5.08e+00 5.13e-01 -3.8 1.01e+02 -4.1 5.03e-01 1.57e-01h 1
36 1.7546433e+01 5.58e+00 7.41e-01 -3.8 2.84e+02 -4.6 2.23e-01 4.17e-02h 1
37 1.6990149e+01 5.97e+00 8.25e-01 -3.8 3.81e+03 -5.1 9.59e-03 2.67e-03f 1
38 1.6599440e+01 5.72e+00 7.56e-01 -3.8 1.02e+02 -4.7 8.15e-01 8.47e-02h 1
39 1.6358942e+01 5.71e+00 7.33e-01 -3.8 3.13e+02 -5.1 1.42e-01 3.16e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 1.6209837e+01 5.65e+00 7.13e-01 -3.8 4.23e+02 -5.6 1.28e-01 2.76e-02h 1
41 1.6129643e+01 5.53e+00 6.92e-01 -3.8 2.36e+02 - 4.36e-01 2.75e-02h 1
42 1.6057004e+01 5.32e+00 6.31e-01 -3.8 4.11e+01 - 1.00e+00 3.83e-02h 1
43 1.6174920e+01 3.75e+00 2.75e-01 -3.8 8.40e+00 - 8.06e-01 2.94e-01h 1
44 1.6492849e+01 3.42e-01 1.47e-02 -3.8 2.57e+00 - 1.00e+00 9.11e-01h 1
45 1.6526505e+01 1.54e-04 6.83e-05 -3.8 1.65e-01 - 1.00e+00 1.00e+00h 1
46 1.6520275e+01 6.33e-05 1.05e-06 -5.7 9.58e-02 - 1.00e+00 1.00e+00h 1
47 1.6520269e+01 6.60e-11 1.88e-11 -5.7 1.24e-04 - 1.00e+00 1.00e+00h 1
48 1.6520192e+01 9.76e-09 1.62e-10 -8.6 1.19e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 48
(scaled) (unscaled)
Objective...............: 1.6520191521080925e+01 1.6520191521080925e+01
Dual infeasibility......: 1.6218784655920614e-10 1.6218784655920614e-10
Constraint violation....: 9.7603845006233314e-09 9.7603845006233314e-09
Complementarity.........: 2.5852537039547041e-09 2.5852537039547041e-09
Overall NLP error.......: 9.7603845006233314e-09 9.7603845006233314e-09
Number of objective function evaluations = 59
Number of objective gradient evaluations = 49
Number of equality constraint evaluations = 59
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 49
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 48
Total CPU secs in IPOPT (w/o function evaluations) = 0.133
Total CPU secs in NLP function evaluations = 0.006
EXIT: Optimal Solution Found.
final time = 16.52
3.4.1.4.3. Plot Results#
# Define empty lists
x = [] # position, units of length
v = [] # velocity, units of length per time
u = [] # acceleration, units of length per time squared
time=[] # time
# Loop over time and append the solution values for each variable to their respective lists
for i in m.tau:
time.append(pyo.value(m.time[i]))
x.append(pyo.value(m.x[i]))
v.append(pyo.value(m.v[i]))
u.append(pyo.value(m.u[i]))
# Make a figure
plt.figure(figsize=(12,4))
# Format subplot 1 (position)
plt.subplot(131)
plt.plot(time,x,linewidth=3,label='x')
plt.title('location (m)',fontsize=16, fontweight='bold')
plt.xlabel('time (s)',fontsize=16, fontweight='bold')
plt.tick_params(direction="in",labelsize=15)
# Format subplot 2 (velocity)
plt.subplot(132)
plt.plot(time,v,linewidth=3,label='v')
plt.xlabel('time (s)',fontsize=16, fontweight='bold')
plt.title('velocity (m/s)',fontsize=16, fontweight='bold')
plt.tick_params(direction="in",labelsize=15)
# Format subplot 3 (acceleration)
plt.subplot(133)
plt.plot(time,u,linewidth=3,label='a')
plt.xlabel('time (s)',fontsize=16, fontweight='bold')
plt.title('acceleration (m/s2)',fontsize=16, fontweight='bold')
plt.tick_params(direction="in",labelsize=15)
plt.show()
Discussion Questions
How does the time/number of evaluations that IPOPT needs to solve the problem change for a discretized versus non-discretized model?
Why might the two solutions differ?
Check to make sure that these results make sense. Based on the known derivative relationships between position, velocity, and acceleration, do the 3 plots make sense? Hint: Match up the trends in each graph.