Tutorial
Formulating an optimal control problem in DYNA is straightforward, and this guide will assist you in setting up the problem correctly and efficiently.
A general purpose optimal control problem can be formulated in the following form \[ \min_{\substack{x,u,p,t_0,t_f}} \Phi\big(x(t_0),t_0,x(t_f),t_f,p\big) + \int_{t_0}^{t_f} L\big(x(t),u(t),t,p\big) \ \mathrm{d}t \] (named Bolza form) subject to \[ \dot{x}(t) = f\big(x(t),u(t),t,p\big), \ \forall \ t \in [t_0,t_f]\\ c_{\text{min}}(t) \leq c\big(x(t),u(t),t,p\big) \ \leq c_{\text{max}}(t), \ \forall \ t \in [t_0,t_f]\\ b_{\text{min}} \leq b\big(x(t_0),t_0,x(t_f),t_f,p\big) \leq b_{\text{max}}\\ p_{\text{min}} \leq p \leq p_{\text{max}} \] with
- \(x(t) \in \mathbb{R}^n\) the state of the system,
- \(u(t) \in \mathbb{R}^m\) the control input,
- \(p \in \mathbb{R}^p\) the static parameters,
- \(t_0 \in \mathbb{R}\) and \(t_f \in \mathbb{R}\) the initial and terminal time,
- \(\Phi\) is the Mayer cost functional (\(\Phi: \mathbb{R}^n \times \mathbb{R} \times \mathbb{R}^n \times \mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R} \)),
- \(L\) is the Lagrange cost functional (\(L: \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R} \)),
- \(f\) is the dynamic constraint (\(f: \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R}^n\)),
- \(c\) is the path constraint (\(c: \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R}^c\)) and
- \(b\) is the boundary condition (\(b: \mathbb{R}^n \times \mathbb{R} \times \mathbb{R}^n \times \mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R}^b\))
Table of content
- Installation
- Tutorial A: Closed Loop Control for a Sliding Mass System
- Tutorial B: Car Problem
- Tutorial C: Goddard Rocket
- Tutorial D: Infinite Horizon Optimal Control Problem
- Tutorial E: Minimum Time-to-Climb of a Supersonic Aircraft
- Tutorial F: Methanol to Hydrocarbons Conversion
- Tutorial G: Optimal Control of a Jojo (Multi-stage OCP)
- Tutorial H: Binary Variant of the Van der Pol Oscillator (MIOCP)
- Tutorial I: Heat Equation (PDE)
- Tutorial J: A Two Reach River Pollution Control (Discrete time OCP)
- Tutorial K: Three Body Problem (ODE Simulation)
- Reference Manual
Installation
Unzip the package and copy all the files preserving the folder structure.
Edit the file .\bin\dyna-setinv.cmd and customize the environment variables (path to GAMS, path to GNUPLOT, etc.).
:: --- Begin of customization ----------------------------------------- :: --- Character set chcp 1252 > nul :: --- Optional components (for development) set rexx= set SevenZ="c:\Program Files\7-Zip\7z.exe" set delphi= set fpc= :: --- Recommended components set gnuplot="c:\wgnuplot506\bin\wgnuplot.exe" set octave="c:\octave-5.1.0.0\mingw64\bin\octave-cli.exe" :: --- Mandatory components set gamsp=C:\GAMS\win64\30.2 set dyna2gams_work=D:\DYNA2GAMS set dyna2gams_home=D:\DYNA2GAMS set dyna-to-gams="%dyna2gams_home%\bin\dyna-to-gams.exe" set dyna-to-html="%dyna2gams_home%\bin\dyna-to-html.exe" set dyna-csv-echelon="%dyna2gams_home%\bin\dyna-csv-echelon.exe" set dyna-csv-reshape="%dyna2gams_home%\bin\dyna-csv-reshape.exe" set dyna-csv-plot="%dyna2gams_home%\bin\dyna-csv-plot.exe" set dyna-csv-splot="%dyna2gams_home%\bin\dyna-csv-splot.exe" :: --- End of customization -------------------------------------------
Execute either .\bin\dyna-setinv.cmd from the Windows command-line environment or .\bin\dyna-shell.cmd from Windows File Explorer. In the console windows, key in dyna-check-install to check everything is fine.
Congratulations ! DYNA2GAMS is now installed.
In all the examples here after, we assume that D:\DYNA2GAMS is the home folder for DYNA2GAMS and the work folder as well.
In the work folder, you can also have a dyna-config.ini file to customize DYNA2GAMS. This is explained in the reference manual.
Tutorial A: Closed Loop Control for a Sliding Mass System
Find \(u\) over \(t \in [0,t_f]\) to minimize \[ J = \int_0^{t_f} (y_1^2 + y_2^2 + u^2) \ \mathrm{d}t \] subject to: \[ \dot{y}_1 = y_2 \\ \dot{y}_2 = u - y_2 |y_2| \\ y_1(0) = 2, \ \ y_2(0) = -2 \] The DYNA implementation is as follows:
par: # declare parameters: here the final time
tf = 1
dyn: # declare the dynamic variables
y{1:2}
J u
t=t0: # set the initial conditions
y{1:2} = {2 -2}
J = 0
equ: # define the equations
y1´ == y2 # character ´ (hexadecimal B4, decimal 180) denotes the derivative
y2´ == u - y2*y2*sign(y2)
J´ == 0.5 * (sqr(y1) + sqr(y2) + sqr(u))
obj: # set the objective and select the solver to use
# (in this case a discontinuous NLP solver)
minimize final(J) using dnlp
How to run this model ? First check that the model syntax is fine, both for DYNA and GAMS.
[Dyna2Gams] D:\DYNA2GAMS>dyna-compile .\examples\tutor-a.dyn Dyna-to-Gams Version 1.0 - Copyright (c) 2018 Alain J. Michiels. All rights reserved. Translating '.\examples\tutor-a.dyn'... Model name: tutor_a Model type: DNLP Number of parameters: 1 Number of dynamic variables: 4 Number of equations: 3 List of parameters: tf List of dynamic variables: y1 y2 J u Calling Gams... *** Status: Normal completion [Dyna2Gams] D:\DYNA2GAMS>dir /w *tutor-a* | find "." $tutor-a.gms $tutor-a.log $tutor-a.lst
The “compilation” checks that the model syntax is fine, both for DYNA and GAMS. It should end by “Normal completion”.
Three files have been created: a .gms file with the model translated in the GAMS language as well as a .log and a .lst created by GAMS during the compilation phase. The compilation step is optional: the sole purpose is to check that the syntax is correct.
The real job is carried out by dyna-execute. By default, four steps are accomplished with a progressive refinement of the mesh grid: 48, 96, 192 and 384 time steps.
[Dyna2Gams] D:\DYNA2GAMS>dyna-execute .\examples\tutor-a.dyn Dyna-to-Gams Version 1.0 - Copyright (c) 2018 Alain J. Michiels. All rights reserved. Translating '.\examples\tutor-a.dyn'... Model name: tutor_a Model type: DNLP Number of parameters: 1 Number of dynamic variables: 4 Number of equations: 3 List of parameters: tf List of dynamic variables: y1 y2 J u Calling Gams... Optimization launched (N=48|CM=L3A4|CF=IL) **** SOLVER STATUS 1 Normal Completion **** MODEL STATUS 2 Locally Optimal **** OBJECTIVE VALUE 1.6173 Optimization launched (N=96|CM=L3A4|CF=IL) **** SOLVER STATUS 1 Normal Completion **** MODEL STATUS 2 Locally Optimal **** OBJECTIVE VALUE 1.6173 Optimization launched (N=192|CM=L3A4|CF=IL) **** SOLVER STATUS 1 Normal Completion **** MODEL STATUS 2 Locally Optimal **** OBJECTIVE VALUE 1.6173 Optimization launched (N=384|CM=L3A4|CF=IL) **** SOLVER STATUS 1 Normal Completion **** MODEL STATUS 2 Locally Optimal **** OBJECTIVE VALUE 1.6173
Let’s have a look at the results. Many files are produced, the user familiar with GAMS will recognize the file types. The one of uttermost interest is the .csv that holds the results.
[Dyna2Gams] D:\DYNA2GAMS>dir /w *tutor-a* | find "." $$tutor-a.192.csv $$tutor-a.192.gdx $$tutor-a.192.log $$tutor-a.192.lst $$tutor-a.192.nfo $$tutor-a.384.csv $$tutor-a.384.gdx $$tutor-a.384.log $$tutor-a.384.lst $$tutor-a.384.nfo $$tutor-a.48.csv $$tutor-a.48.log $$tutor-a.48.lst $$tutor-a.48.nfo $$tutor-a.768.gdx $$tutor-a.96.csv $$tutor-a.96.gdx $$tutor-a.96.log $$tutor-a.96.lst $$tutor-a.96.nfo $tutor-a.gms $tutor-a.log $tutor-a.lst [Dyna2Gams] D:\DYNA2GAMS>%gamsp%\gbin\head $$tutor-a.48.csv "_N_","Time","y1","y2","J","u" "N0", 0.000000E+00, 2.000000E+00,-2.000000E+00, 0.000000E+00, 1.548745E-01 "N1", 1.727458E-02, 1.966057E+00,-1.930623E+00, 6.752117E-02, 1.547967E-01 "N2", 4.522542E-02, 1.913558E+00,-1.827664E+00, 1.697485E-01, 1.550753E-01 "N3", 6.250000E-02, 1.882497E+00,-1.769126E+00, 2.289988E-01, 1.554874E-01 "N4", 7.977458E-02, 1.852417E+00,-1.714051E+00, 2.855201E-01, 1.560419E-01 "N5", 1.077254E-01, 1.805679E+00,-1.631507E+00, 3.716928E-01, 1.571662E-01 "N6", 1.250000E-01, 1.777908E+00,-1.584138E+00, 4.219597E-01, 1.579831E-01 "N7", 1.422746E-01, 1.750934E+00,-1.539277E+00, 4.701264E-01, 1.588597E-01 "N8", 1.702254E-01, 1.708870E+00,-1.471503E+00, 5.439556E-01, 1.603422E-01
Tutorial B: Car Problem
The classical optimal control example, the car problem, is considered. The objective is to find an optimal acceleration strategy for a car to travel a predetermined distance in minimum time. The car starts and stops at zero velocity and it is assumed that there is no traveling speed limit.
The problem is to find \(u\) over \(t \in [0,t_f]\) that minimizes \[ J = t_f \] subject to: \[ \dot{x}_1 = u \\ \dot{x}_2 = x_1 \\ -2 \leq u \leq 1 \\ x_1(0) = 0, \ \ x_2(0) = 0 \\ x_1(t_f) = 0, \ \ x_2(t_f) = 300 \] The model can be solved analytically and gives the optimal solution \(t_f^* = 30\).
The DYNA implementation reads as follows:
var: # declare the time horizon as a static variable tf dyn: # declare the dynamic variables x1 x2 u lim: # define lower and upper limits -2 <= u <= 1 0 <= tf <= 50 t=t0: # set the initial conditions x1 = 0 x2 = 0 t=tf: # set the final conditions x1 = 0 x2 = 300 equ: # define the equations of motion x1´ == u x2´ == x1 obj: # set the objective and select the solver to # use (in this case the non linear solver IPOPT) minimize tf using nlp with ipopt
A reasonable upper limit is set on \(t_f\) for practical reasons. We could define \(t_f\) on the interval \([0, +\infty)\), in which case we need to provide the solver with a somewhat appropriate initial guess of its optimal value. By default, DYNA uses the bounds on the variables to generate initial guesses. If an upper and a lower bound is given, the initial guess will be the arithmetic mean. If only one of these bounds is specified (while the other one is \(\pm \infty\)) the initial guess will be equal to this bound. If there is a variable for which no bounds are specified, the initial guess will simply be 0.
var: # declare the time horizon as a static variable tf dyn: # declare the dynamic variables x1 x2 u lim: # define lower and upper limits -2 <= u <= 1 0 <= tf <= +inf t=t0: # set the initial conditions x1 = 0 x2 = 0 t=tf: # set the final conditions x1 = 0 x2 = 300 ini: # guess an appropriate starting point tf = 25 equ: # define the equations of motion x1´ == u x2´ == x1 obj: # set the objective and select the solver to # use (in this case the non linear solver IPOPT) minimize tf using nlp with ipopt
Tutorial C: Goddard Rocket
Maximize the final altitude of a vertically launched rocket, using the thrust as a control and given the initial mass, the fuel mass, and the drag characteristics of the rocket.
Find \(T\) over \(t \in [0,t_f]\) to maximize \[ J = h(t_f) \] subject to: \[ \dot{h} = v \\ m \ \dot{v} = T - D(v,h) - m \ g(h) \\ \dot{m} = \frac{-T}{c} \\ D(v,h) = D_c v^2 e^{-h_c \big(\frac{h-h_0}{h_0}\big)} \\ g(h) = g_0 \big(\frac{h_0}{h}\big)^2 \\ h(0) = h_0, \ \ v(0) = v_0, \ \ m(0) = m_0 \\ m(t_f) = m_f \\ 0 \leq T \leq T_{\text{max}} \]
par: # define the parameters of the model v0 = 0.0 m0 = 1.0 g0 = 1.0 h0 = 1.0 hc = 500.0 vc = 620.0 mc = 0.6 Tmax = 3.5*g0*m0 c = 0.5*sqrt(g0*h0) mf = mc*m0 Dc = 0.5*vc*m0/g0 var: # the time horizon is an endogenous variable tf dyn: # declare the dynamic variables h v m T lim: # box constraint the variables 0.01 <= tf <= 1 h0 <= h <= +inf v0 <= v <= +inf mf <= m <= m0 0 <= T <= Tmax t=t0: # set the initial conditions h = h0 v = v0 m = m0 ini: # provide the solver with own initial guesses tf = 1.0 h = h0 v = v0 m = m0 T = Tmax/2 exp: # define the symbolic variables D == Dc*sqr(v)*exp(-hc*(h-h0)/h0) g == g0*sqr(h0/h) equ: # define the equations h´ == v v´ == (T - D)/m - g m´ == -T/c obj: # set the objective and select the solver to # use (in this case a continuous NLP solver) maximize final(h) using nlp
Tutorial D: Infinite Horizon Optimal Control Problem
Consider the following optimal control problem taken from [1]. Denoting \(x(t) = [x_1(t) \ x_2(t)]^T \in \mathbb{R}^2\) as the state and \(u(t) \in \mathbb{R}\) as the control, minimize the cost functional \[ J = \frac{1}{2} \int_0^{\infty} (x^TQx + u^TRu) \ \mathrm{d}t \] subject to the dynamic constraint \[ \dot x = Ax + Bu \] and the initial condition \[ x(0) = \begin{bmatrix} -4 \\ 4 \end{bmatrix} \] The matrices \(A, B, Q,\) and \(R\) for this problem are given as \[ A = \begin{bmatrix} 0 & 1 \\ 2 & -1 \end{bmatrix} \quad,\quad B = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \quad,\quad Q = \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix} \quad,\quad R = \frac{1}{2} \] The exact solution to this problem is \[ x(t) = \exp([A-BK]t) \ x(0) \\ u(t) = -Kx(t) \\ \lambda(t) = Sx(t) \] where \(K\) is the optimal feedback gain and \(S\) is the solution to the algebraic Riccati equation. In this case \(K\) and \(S\) are given, respectively, as \[ K = \begin{bmatrix} 4.8284 & 2.5576 \end{bmatrix} \\ S = \begin{bmatrix} 6.0313 & 2.4142 \\ 2.4142 & 1.2788 \end{bmatrix} \]
References:
- F. Fahroo and I.M. Ross, Pseudospectral Methods for Infinite-Horizon Nonlinear Optimal Control Problems, Journal of Guidance, Control, and Dynamics, 31(4):927-936, 2008.
- D. Garg, W.W. Hager, A.V. Rao, Gauss Pseudospectral Method for Solving Infinite-Horizon Optimal Control Problems, AIAA Guidance, Navigation, and Control Conference 2-5 August 2010, Toronto, Ontario Canada
par: # define the parameters
A{11 12 21 22} = {0 1 2 -1}
B{1 2} = {0 1}
Q{11 12 21 22} = {2 0 0 1}
R = 0.5
K{1 2} = {4.8284 2.5576}
S{11 12 21 22} = {6.0313 2.4142 2.4142 1.2788}
t{1:2} = {0.83378 0.42404}
Y{11 12 21 22} = {-4.68805 0.68805 5.62263 -1.62263}
tf = 99
var: # define the static variables
uK{1:2}
dyn: # declare the dynamic variables
t
x{1:2}
u
J
lambda{1:2}
y{1:2}
v
mu{1:2}
lim: # set lower and upper bounds
1 <= uK{1:2} <= 10
0 <= t <= +inf
t=t0: # set the initial conditions
t = 0
x{1:2} = {-4 4}
J = 0
t=tf: # set the final conditions
lambda{1:2} = 0
equ: # define the equations of the system
t´ == 1
# numerical solution follows
e1.. x1´ == (A11*x1 + A12*x2 + B1*u)
e2.. x2´ == (A21*x1 + A22*x2 + B2*u)
J´ == 0.5*(Q11*sqr(x1) + (Q12+Q21)*x1*x2 + Q22*sqr(x2) + R*sqr(u))
u == - uK1*x1 - uK2*x2
lambda1´ == - Q11*x1 - Q12*x2 - A11*lambda1 - A21*lambda2
lambda2´ == - Q21*x1 - Q22*x2 - A12*lambda1 - A22*lambda2
# analytical solution follows
y1 == Y11*exp(-t/t1) + Y12*exp(-t/t2)
y2 == Y21*exp(-t/t1) + Y22*exp(-t/t2)
v == - K1*y1 - K2*y2
mu1 == S11*y1 + S12*y2
mu2 == S21*y1 + S22*y2
obj: # set the objective and select the solver to use
minimize final(J) using nlp
gms: # a bit of GAMS code that will be generated by DYNA
@costates # calculate the costates from the duals
@tpa eta1 = costate(e1) # retrieve costate 1
@tpa eta2 = costate(e2) # retrieve costate 2
@csvsave eta1 eta2 # save them in the results file
gpl: # visualize the results with GNUPLOT
set multiplot layout 2,2 rowsfirst
@plotYYT 2, x1, y1, "States 1 vs. Time"
@plotYYT 2, x2, y2, "States 2 vs. Time"
@plotYYT 3, eta1, lambda1, mu1, "Costates 1 vs. Time"
@plotYYT 3, eta2, lambda2, mu2, "Costates 2 vs. Time"
unset multiplot
Now, using GNUPLOT...
[Dyna2Gams] D:\DYNA2GAMS>dyna-gnuplot tutor-d 384 Invoking Gnuplot... Press any key to continue...
Tutorial E: Minimum Time-to-Climb of a Supersonic Aircraft
In this example, we will introduce the tabulated functions. The main purpose of tabulated functions is to represent functions that don’t have an algebraic form. DYNA will perform an interpolation using various schemes ranging from linear interpolation to cubic splines, preserving shapes or not, smooth or non-smooth.
The minimum time-to-climb of a supersonic aircraft optimal control problem is stated as follows.
Minimize the cost functional \[ J = t_f \] subject to the dynamic constraints \[ \dot{h}(t) = v(t)\sin\alpha(t) \\ \dot{v}(t) = \frac{T\cos\alpha(t)-D}{m(t)} \\ \dot{\gamma}(t) = \frac{T\sin\alpha(t)+L}{m(t)v(t)} + \big(\frac{v(t)}{r(t)} - \frac{\mu}{v(t)r^2(t)}\big) \cos\gamma(t) \\ \dot{m}(t) = -\frac{T}{g_0 I_{sp}} \] and the boundary conditions \[ h(0)=0, \ h(t_f)=19994.88 \\ v(0)=129.314, \ v(t_f)=295.092 \\ \gamma(0)=0, \ \gamma(t_f)=0 \\ m(0)=19050.864 \] This example is taken verbatim from the following reference:
- A.E. Bryson, M.N. Desai, W.C. Hoffman, Energy-State Approximation in Performance Optimization of Supersonic Aircraft, Journal of Aircraft, 6(6):481-488, November-December 1969.
fun: # character § (hexadecimal A7, decimal 167) denotes a DYNA tabulated function
# Bryson Minimum Time To Climb Aero Data
# U.S. 1976 Standard Atmosphere
# Column 1: Altitude (m)
# Column 2: Atmospheric Density (kg/m³)
# Column 3: Speed of Sound (m/s)
ALT §AD §SS
-2000 1.478e+00 3.479e+02
0 1.225e+00 3.403e+02
2000 1.007e+00 3.325e+02
4000 8.193e-01 3.246e+02
6000 6.601e-01 3.165e+02
8000 5.258e-01 3.081e+02
10000 4.135e-01 2.995e+02
12000 3.119e-01 2.951e+02
14000 2.279e-01 2.951e+02
16000 1.665e-01 2.951e+02
18000 1.216e-01 2.951e+02
20000 8.891e-02 2.951e+02
22000 6.451e-02 2.964e+02
24000 4.694e-02 2.977e+02
26000 3.426e-02 2.991e+02
28000 2.508e-02 3.004e+02
30000 1.841e-02 3.017e+02
32000 1.355e-02 3.030e+02
34000 9.887e-03 3.065e+02
36000 7.257e-03 3.101e+02
38000 5.366e-03 3.137e+02
40000 3.995e-03 3.172e+02
42000 2.995e-03 3.207e+02
44000 2.259e-03 3.241e+02
46000 1.714e-03 3.275e+02
48000 1.317e-03 3.298e+02
50000 1.027e-03 3.298e+02
52000 8.055e-04 3.288e+02
54000 6.389e-04 3.254e+02
56000 5.044e-04 3.220e+02
58000 3.962e-04 3.186e+02
60000 3.096e-04 3.151e+02
62000 2.407e-04 3.115e+02
64000 1.860e-04 3.080e+02
66000 1.429e-04 3.044e+02
68000 1.091e-04 3.007e+02
70000 8.281e-05 2.971e+02
72000 6.236e-05 2.934e+02
74000 4.637e-05 2.907e+02
76000 3.430e-05 2.880e+02
78000 2.523e-05 2.853e+02
80000 1.845e-05 2.825e+02
82000 1.341e-05 2.797e+02
84000 9.690e-06 2.769e+02
86000 6.955e-06 2.741e+02
# Propulsion Data for Bryson Aircraft
# The thrust for the aircraft considered by Bryson in 1969 depends
# upon the Mach number and the altitude. This data is taken verbatim
# from the 1969 Journal of Aircraft paper (see reference above) and
# is copied for use in this example. The data are stored in the
# following table:
# Column 1: Mach number
# Row 1 : altitude (meters)
# Data : Thrust (Newton)
§Thrust = GridXY
0 1524 3048 4572 6096 7620 9144 12192 15240 21336
0 107646.9724 106757.328 90298.9066 76954.2406 64499.219 54268.3084 45371.8644 25354.8654 15123.9548 444.8222
0.2 124550.216 109426.2612 93857.4842 80512.8182 67612.9744 56937.2416 47595.9754 28913.443 17348.0658 889.6444
0.4 125884.6826 112095.1944 97416.0618 83181.7514 70726.7298 59606.1748 49820.0864 32472.0206 19572.1768 1779.2888
0.6 137005.2376 120991.6384 105867.6836 91188.551 76954.2406 65388.8634 54713.1306 36030.5982 21796.2878 3558.5776
0.8 153463.659 134781.1266 118322.7052 103198.7504 88074.7956 74730.1296 62719.9302 41813.2868 24910.0432 4893.0442
1 168587.6138 152574.0146 135225.9488 119212.3496 103643.5726 88074.7956 74730.1296 49820.0864 30247.9096 6227.5108
1.2 160580.8142 169032.436 155242.9478 139229.3486 121436.4606 104978.0392 89409.2622 59606.1748 36920.2426 7561.9774
1.4 160580.8142 162804.9252 171256.547 160580.8142 140563.8152 124995.0382 107646.9724 72061.1964 44482.22 9786.0884
1.6 160580.8142 156577.4144 187270.1462 172146.1914 158801.5254 142343.104 124995.0382 85850.6846 52933.8418 12899.8438
1.8 160580.8142 150349.9036 203283.7454 183711.5686 177039.2356 153908.4812 138339.7042 96526.4174 59161.3526 13789.4882
# Aerodynamic Data for Bryson Aircraft
# M: Mach number values
# CLalpha: coefficient of lift values
# CD0: zero-lift coefficient of drag values
# eta: load factors
M §CLalpha §CD0 §eta
0.0 3.44 0.013 0.54
0.4 3.44 0.013 0.54
0.8 3.44 0.013 0.54
0.9 3.58 0.014 0.75
1.0 4.44 0.031 0.79
1.2 3.44 0.041 0.78
1.4 3.01 0.039 0.89
1.6 2.86 0.036 0.93
1.8 2.44 0.035 0.93
par:
Re = 6378145
mu = 3.986e14
S = 49.2386
g0 = 9.80665
Isp = 1600
var:
tf
dyn:
h v fpa m # state variables
alpha # control variable
lim:
100 <= tf <= 600
0 <= h <= 21031.2
5 <= v <= 1000
-40/180*pi <= fpa <= +40/180*pi
22 <= m <= 20410
-45/180*pi <= alpha <= +45/180*pi
t=t0:
h = 0 # Initial altitude (meters)
v = 129.314 # Initial speed (m/s)
fpa = 0 # Initial flight path angle (rad)
m = 19050.864 # Initial mass (kg)
t=tf:
h = 19994.88 # Final altitude (meters)
v = 295.092 # Final speed (m/s)
fpa = 0 # Final flight path angle (rad)
exp:
CD == CD0 + eta * CLalpha * alpha * alpha
CL == CLalpha * alpha
q == 0.5 * rho * v * v
D == q * S * CD # D is the magnitude of the drag force
L == q * S * CL # L is the magnitude of the lift force
r == h + Re
rho == §AD(h)
sos == §SS(h)
Mach == v / sos
CD0 == §CD0(Mach)
CLalpha == §CLalpha(Mach)
eta == §eta(Mach)
Thrust == §Thrust(Mach,h)
ini:
tf = 300
h = linspace(initial(h),final(h))
v = linspace(initial(v),final(v))
fpa = linspace(initial(fpa),final(fpa))
m = initial(m)
alpha = linspace(45/180*pi,-45/180*pi)
equ:
h´ == v * sin(fpa)
v´ == (Thrust * cos(alpha) - D) / m - mu * sin(fpa) / (r * r)
fpa´ == (Thrust * sin(alpha) + L) / (m * v) + cos(fpa) * (v / r - mu / (v * r * r))
m´ == -Thrust / (g0 * Isp)
obj:
minimize tf using nlp
Tutorial F: Methanol to Hydrocarbons Conversion
In this example, we introduce tabulated observations for a dynamic system to identify its parameters.
Determine the reaction coefficients for the conversion of methanol into various hydrocarbons.
The nonlinear model that describes the process is \[ \dot{x}_1 = - \big(2\theta_2 - \frac{\theta_1 x_2}{(\theta_2+\theta_5)x_1+x_2} + \theta_3 + \theta_4 \big) x_1 \\ \dot{x}_2 = \frac{\theta_1 x_1(\theta_2 x_1-x_2)}{(\theta_2+\theta_5)x_1+x_2} + \theta_3 x_1 \\ \dot{x}_3 = \frac{\theta_1 x_1(x_2+\theta_5 x_1)}{(\theta_2+\theta_5)x_1+x_2} + \theta_4 x_1 \\ \] with coefficients \(\theta_i \geq 0\). Initial conditions are known.
See: COPS: Large-Scale Optimization Problems
Reference:
- E.D. Dolan, J.J. Moré and T.S. Munson, Benchmarking Optimization Software with COPS 3.0, Argonne National Laboratory, Technical Report ANL/MCS-TM-273, February 2004.
par:
tf = 1.122
var:
theta{1:5}
dyn:
x{1:3}
obs: # tabulation of the observations
Time x1 x2 x3
0.000 1.0000 0.0000 0.0000
0.050 0.7085 0.1621 0.0811
0.065 0.5971 0.1855 0.0965
0.080 0.5537 0.1989 0.1198
0.123 0.3684 0.2845 0.1535
0.233 0.1712 0.3491 0.2097
0.273 0.1198 0.3098 0.2628
0.354 0.0747 0.3576 0.2467
0.397 0.0529 0.3347 0.2884
0.418 0.0415 0.3388 0.2757
0.502 0.0261 0.3557 0.3167
0.553 0.0208 0.3483 0.2954
0.681 0.0085 0.3836 0.2950
0.750 0.0053 0.3611 0.2937
0.916 0.0019 0.3609 0.2831
0.937 0.0018 0.3485 0.2846
1.122 0.0006 0.3698 0.2899
lim:
0 <= theta{1:5} <= 3
0 <= x{1:3} <= 2
t=t0:
x1 = 1.0
x{2:3} = 0.0
ini:
x{1:3} = ObservedVal(x{1:3}) # a dedicated function to get access to the data
theta{1:5} = 1
equ:
x1´ == -(2*theta2-(theta1*x2)/((theta2+theta5)*x1+x2)+theta3+theta4)*x1
x2´ == (theta1*x1*(theta2*x1-x2))/((theta2+theta5)*x1+x2)+theta3*x1
x3´ == (theta1*x1*(x2+theta5*x1))/((theta2+theta5)*x1+x2)+theta4*x1
obj: # a dedicated function for the lack-of-fit
minimize sum{1:3, ObservedSSE(x{:})} using nlp
gms: # dedicated keywords to build a customized report
parameter report(*,*);
set o(Observation);
ObservedSet(x{1:3},o);
$$ondotL
report(o,'x{1:3}') = ObservedGap(x{1:3},o);
display report;
Tutorial G: Optimal Control of a Jojo (Multi-stage OCP)
As an example for a multi-stage optimization problem, we consider the simple jojo described in [1]. Here, the altitude of the jojo is denoted by the state variable \(y\), while \(w = dy/dt\) is the associated velocity. The hand playing with the jojo is described by the altitude \(x\), while \(v = dx/dt\) is the velocity of the hand. The distance \(l\) between hand and jojo is defined as \(l = x - y\). As soon as the jojo is released from the hand, the equations of motion are given by \[ \dot{x}(t) = v(t) \\ \dot{v}(t) = u(t) \\ \dot{y}(t) = w(t) \\ \dot{w}(t) = -\mu g + k u + a(w-v) \] Here, we define the jojo’s damping ratio \(k = J / (m r^2 + J)\) as well as the effective mass ratio \(\mu = 1 - k\), which depend on the inertia \(J\), the mass \(m\) and the (inner) coiling radius \(r\) of the jojo. Moreover, \(g\) denotes the acceleration due to gravity while \(a\) is the jojo’s coiling friction. Note that we assume here that the hand can directly be controlled by the input function \(u\).
Let us assume that the sum of \(r\) and the total length of the rope is given by \(L\). When the distance \(l = x - y\) happens to be equal to this length \(L\), the kinetic energy of the jojo is suddenly absorbed, while for the rotational energy a conservation law is satisfied. Thus, the jump is given by \[ x(t^+_j) = x(t^-_j) \\ v(t^+_j) = v(t^-_j) \\ y(t^+_j) = y(t^-_j) \\ w(t^+_j) = k v(t^-_j) + \sqrt{k^2 v^2(t^-_j) + k(w^2(t^-_j) - 2 w(t^-_j) v(t^-_j))} \] In order to understand this expression, it is helpful to remark that for the passive case \(v = 0\) (i.e. for the case that the hand is at rest during the jump) the jump condition simplifies to \[ w(t^+_j) = \sqrt{k} \ w(t^-_j) \] i.e. the damping ratio \(k\) measures the loss of kinetic energy for a passive jump. After the jump, the equations of motion are again given by the above dynamic system. Now, our aim is to minimize the control action of the hand defined as \[ \Phi = \int_0^{t_f} u^2(t) \ \textrm{d}t \] while starting and stopping the jojo under the following conditions: \[ x(0) = y(0) = 0 \ \textrm{m}\\ v(0) = w(0) = 0\ \textrm{m/s} \\ x(t_j) - y(t_j) = L \\ x(t_f) = y(t_f) = -0.1\ \textrm{m} \\ v(t_f) = w(t_f) = 0\ \textrm{m/s} \] i.e. we start both the hand and the jojo at rest at the position 0. Moreover at the time \(t_j\) we require that the length of the rope is equal to \(L = 1\) m. Finally, the jojo should be stopped when softly touching the hand \(-10\) cm below the starting point. Note that the durations of the two phases, represented by the time variables \(t_j\) and \(t_f-t_j\), are optimized, too.
Note that the velocity of the hand has a peek at the time, where the jump is. The velocity of the jojo is discontinuous at this time. One interpretation of this result is that the energy loss during the jump can be reduced by moving the hand quickly upwards during this phase.
Reference:
- B. Houska, H.J. Ferreau and M. Diehl, ACADO Toolkit - An open-source framework for automatic control and dynamic optimization, Optimal Control Applications and Methods, 32(3):298-312, 2011.
set: phase = ph1:ph2 par: m = 0.200 # the mass of the jojo [kg] J = 1e-4 # the inertia of the jojo [kg*m²] r = 0.010 # the coiling radius of the jojo [m] g = 9.81 # the gravitational constant [m/s²] a = 1e-2 # the coiling friction [1/s] L = 1.00 # the length of the rope [m] k = J/(m*r*r+J) # the jojo’s damping ratio mu = 1.0 - k # the effective mass ratio tf = 1 # simulation time var: tlen[phase] # the duration of the phases dyn: t[phase] # the real time x[phase] # the position of the “hand” v[phase] # the velocity of the “hand” y[phase] # the position of the jojo w[phase] # the velocity of the jojo u[phase] :: control # the control action of the “hand” phi[phase] # the cost functional lim: 0.1 <= tlen[phase] <= 2 -10 <= u[phase] <= 10 t=t0: t['ph1'] = 0 x['ph1'] = 0 v['ph1'] = 0 y['ph1'] = 0 w['ph1'] = 0 phi['ph1'] = 0 t=tf: x['ph2'] = -0.1 v['ph2'] = 0 y['ph2'] = -0.1 w['ph2'] = 0 exp: TimeDot == tlen[phase] / tf mac: @ssqrt(x) sqrt(sqrt(sqr(&x))) equ: t´[phase] / TimeDot == 1 x´[phase] / TimeDot == v[phase] v´[phase] / TimeDot == u[phase] y´[phase] / TimeDot == w[phase] w´[phase] / TimeDot == -mu*g + k*u[phase] + a*(v[phase]-w[phase]) phi´[phase] / TimeDot == sqr(u[phase]) @linkphases -x w # link final and initial conditions across phases for state and control variables, except for variable w, which is given here after [phase]$(ord(phase)<card(phase)).. 0 == - initial(w[phase+1]) + final(k*v[phase] + @ssqrt(sqr(k*v[phase]) + k*(sqr(w[phase])-2*w[phase]*v[phase]))) L == final(x['ph1']) - final(y['ph1']) obj: minimize final(phi['ph2']) using nlp
Postprocess the csv file so as to stack the results vertically.
[Dyna2Gams] D:\DYNA2GAMS>dyna-csv-phases tutor-g 384 -2 Dyna-CSV-Reshape Version 1.0 - Copyright (c) 2018 Alain J. Michiels. All rights reserved. Reshaping '$$tutor-g.384.csv'... Header: "_N_","Time","t[ph1] ... phi[ph1]","phi[ph2]" Phase 1: "_N_","t[ph1]","x[ph ... ,"u[ph1]","phi[ph1]" Phase 2: "_N_","t[ph2]","x[ph ... ,"u[ph2]","phi[ph2]" File '$$tutor-g.384.rsh.csv' created.
Tutorial H: Binary Variant of the Van der Pol Oscillator (MIOCP)
This case describes a Van der Pol Oscillator variant with three binary controls instead of only one continuous control.
The mixed-integer optimal control problem is given by \[ \min\limits_{x,y,w} \int_{t_0}^{t_f} \big(x^2(t)+y^2(t)\big) \ \mathrm{d}t \\ \] subject to \[ \dot x = y \\ \dot y = \sum\limits_{i=1}^{3} c_{i} \; w_i \; (1-x^2) y-x \\ x(t_0) = 1 \\ y(t_0) = 0 \\ 1 = \sum\limits_{i=1}^{3}w_i(t) \\ w_i(t) \in \{0,1\} \quad i=1 \ldots 3 \] Source: mintoc.de
Reference:
- S. Sager, M. Jung, C. Kirches, Combinatorial integral approximation, Mathematical Methods of Operations Research, 73:363, June 2011.
set: i = i1:i3 par: c[i] = |i1 -1, i2 0.75, i3 -2| tf = 20 dyn: x y w[i] J lim: 0 <= w[i] <= 1 t=t0: x = 1 y = 0 J = 0 ini: x = 1 y = 0 equ: x´ == y y´ == sum(i, c[i]*w[i])*(1-x*x)*y - x J´ == x*x + y*y $all.. 1 == sum(i, w[i]) # enforce the constraint on all time steps of the horizon obj: minimize final(J) using nlp gms: # Combinatorial Integral Approximation $$ifthen %ITER%==%ITERMAX% @seekcia ControlVar=w[i], Optcr=0.5, Method=plain $$endif
Tutorial I: Heat Equation (PDE)
The heat equation is a parabolic partial differential equation that describes the distribution of heat (or variation in temperature) in a given region over time.
In this example, we seek \(z(x,t)\) over \(x \in [0,1]\) and \(t \in [0,1]\) such that \[ \frac{\partial{z}}{\partial{t}} = D \ \nabla^2 z \\ z(0,t) = 0 \\ z(1,t) = 0 \\ z(x,0) = a \sin \pi x \]
set: x = |x0:x20| xl[x] = |x0| xx[x] = |x1:x19| xr[x] = |x20| par: DeltaX = 1/(card(x)-1) TickX[x] = DeltaX*(ord(x)-1) D = 1 a = 1 tf = 1 dyn: z[x] t=t0: z[x]$xx(x) = a*sin(pi*TickX[x]) equ: $xl[x].. z[x] == 0 $xx[x].. z´[x] == D * (z[x-1]-2*z[x]+z[x+1])/sqr(DeltaX) $xr[x].. z[x] == 0 obj: minimize initial(z['x0']) using lp # dummy objective gpl: @splotZXT z, TickX, "@TITLE@", "z[x,t]", "x"
Visualizing the results
[Dyna2Gams] D:\DYNA2GAMS>dyna-gnuplot tutor-i 384 Dyna-CSV-SPlot Version 1.0 - Copyright (c) 2018 Alain J. Michiels. All rights reserved. Reading '$tutor-i.gpl'... Reading '$$tutor-i.384.gpl'... Reading '$$tutor-i.384.csv'... Header: "_N_","Time","z[x0]" ... ]","z[x19]","z[x20]" File 'gnuplot.csv' created. File 'gnuplot.gpl' created. Invoking Gnuplot... Press any key to continue...
Tutorial J: A Two Reach River Pollution Control (Discrete time OCP)
Simulation of a two reach river pollution control model. The control problem is to maintain the instream biochemical demand (B.O.D.) and the dissolved oxygen (D.O.) at prescribed levels for a river with multiple polluters using the percentage of the B.O.D. treated at the sewage works as the control variable.
Reference:
- M.G. Singh, Dynamical Hierarchical Control, North-Holland, pp. 135-137, 1980.
run:
time-model = Discrete
N = 20
par:
ts = 1
var:
k11 k12 k13 k14 d1
k21 k22 k23 k24 d2
dyn:
x1 x3 # B.O.D. concentration (mg/l) in the stream
x2 x4 # D.O. concentration (mg/l)
u1 u2 # control variables
J
k=k0:
x{1:3} = 0
x4 = 1
J = 0
equ:
# state variables
x1´ == 0.18 * x1 + 2.0 * u1 + 4.50
x2´ == -0.25 * x1 + 0.27 * x2 + 6.15
x3´ == 0.55 * x1 + 0.18 * x3 - 2.0 * u2 + 2.00
x4´ == 0.55 * x2 - 0.25 * x3 + 0.27 * x4 + 2.65
# control law
u1 == k11 * x1 + k12 * x2 + k13 * x3 + k14 * x4 + d1
u2 == k21 * x1 + k22 * x2 + k23 * x3 + k24 * x4 + d2
# criterion
J´ == J + 0.5 * (sqr(x1-5) + sqr(x2-7) + sqr(x3-5) + sqr(x4-7) + 100*sqr(u1) + 100*sqr(u2))
obj:
minimize final(J) using nlp with conopt|snopt
Tutorial K: Three Body Problem (ODE Simulation)
We illustrate the ODE simulation with an example from Astronomy, the restricted three body problem. One considers two bodies of masses \(1-\mu\) and \(\mu\) in circular rotation in a plane and a third body of negligible mass moving around in the same plane. The equations are as follows: \[ \ddot{y}_1 = y_1 + 2 \dot{y}_2 - \mu'\ \frac{y_1+\mu}{D_1} - \mu\ \frac{y_1-\mu'}{D_2} \\ \ddot{y}_2 = y_2 - 2 \dot{y}_1 - \mu'\ \frac{y_2}{D_1} - \mu\ \frac{y_2}{D_2} \\ D_1 = ((y_1+\mu)^2 + y_2^2)^{3/2} \\ D_2 = ((y_1-\mu')^2+y_2^2)^{3/2} \\ \mu = 0.012277471 \\ \mu' =1-\mu \]
Reference:
- E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I - Nonstiff Problems, Springer Series in Computational Mathematics, page 129.
par:
mu = 0.012277471
muc = 1 - mu
tf = 2*17.06521656016
dyn:
y{1:4}
t=t0:
y1 = 0.994
y2 = 0.0
y3 = 0.0
y4 = -2.001585106379
exp:
D1 == (sqr(y1+mu) + sqr(y2))^(3/2)
D2 == (sqr(y1-muc) + sqr(y2))^(3/2)
sim: # Section to describe the ODE problem
y1´ := y3
y2´ := y4
y3´ := y1 + 2*y4 - muc*(y1+mu)/D1 - mu*(y1-muc)/D2
y4´ := y2 - 2*y3 - muc*y2/D1 - mu*y2/D2
Let’s run this model.
[Dyna2Gams] D:\DYNA2GAMS>dyna-execute .\examples\tutor-k-ode.dyn Dyna-to-Gams Version 1.0 - Copyright (c) 2018 Alain J. Michiels. All rights reserved. Translating '.\examples\tutor-k-ode.dyn'... Model name: tutor_k Number of parameters: 3 Number of symbolic variables: 2 Number of dynamic variables: 4 List of parameters: mu muc tf List of symbolic variables: D1 D2 List of dynamic variables: y1 y2 y3 y4 Calling Gams... Simulation launched (N=48|SM=DP853|SF=SH) **** SIMULATION STATUS 1 Solving looks successful **** BENCHMARK VALUE -1.8770 **** MODEL EVALUATIONS 2922 Simulation launched (N=96|SM=DP853|SF=SH) **** SIMULATION STATUS 1 Solving looks successful **** BENCHMARK VALUE 1.9880 **** MODEL EVALUATIONS 2991 Simulation launched (N=192|SM=DP853|SF=SH) **** SIMULATION STATUS 1 Solving looks successful **** BENCHMARK VALUE 0.5990 **** MODEL EVALUATIONS 3564 Simulation launched (N=384|SM=DP853|SF=SH) **** SIMULATION STATUS 1 Solving looks successful **** BENCHMARK VALUE 0.0260 **** MODEL EVALUATIONS 5904 [Dyna2Gams] D:\DYNA2GAMS>dir /w *tutor-k-ode* | find "." $$tutor-k-ode.192.csv $$tutor-k-ode.192.log $$tutor-k-ode.192.lst $$tutor-k-ode.192.nfs $$tutor-k-ode.384.csv $$tutor-k-ode.384.log $$tutor-k-ode.384.lst $$tutor-k-ode.384.nfs $$tutor-k-ode.48.csv $$tutor-k-ode.48.log $$tutor-k-ode.48.lst $$tutor-k-ode.48.nfs $$tutor-k-ode.96.csv $$tutor-k-ode.96.log $$tutor-k-ode.96.lst $$tutor-k-ode.96.nfs $tutor-k-ode.gms $tutor-k-ode.log $tutor-k-ode.lst
The “benchmark value” is the sum of the end values of the state variables. It allows to compare the accuracy of the integration across the simulations. By using finer time grids, we expect that the benchmark values would converge to some number, which is not the case and unfortunately denotes some inaccuracy in the integration.
In the versions tutor-k-dae and tutor-k-dae-multi-phase, we explicitly look for a periodic solution. The simulation of the ODE version (tutor-k-ode) provides the initial guess. The objective of the non linear program is the period of the solution.
Reference Manual
DYNA language comprises a fairly rich set of statements to simplify the description and the solve of optimal control problems. DYNA syntax can be browsed here below along with the description of the DYNA2GAMS suite of tools.