1#
2#     This file is part of CasADi.
3#
4#     CasADi -- A symbolic framework for dynamic optimization.
5#     Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl,
6#                             K.U. Leuven. All rights reserved.
7#     Copyright (C) 2011-2014 Greg Horn
8#
9#     CasADi is free software; you can redistribute it and/or
10#     modify it under the terms of the GNU Lesser General Public
11#     License as published by the Free Software Foundation; either
12#     version 3 of the License, or (at your option) any later version.
13#
14#     CasADi is distributed in the hope that it will be useful,
15#     but WITHOUT ANY WARRANTY; without even the implied warranty of
16#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17#     Lesser General Public License for more details.
18#
19#     You should have received a copy of the GNU Lesser General Public
20#     License along with CasADi; if not, write to the Free Software
21#     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22#
23#
24from casadi import *
25
26"""
27This example mainly intended for CasADi presentations. It contains a compact
28implementation of a direct single shooting method for DAEs using a minimal
29number of CasADi concepts.
30
31It solves the following optimal control problem (OCP) in differential-algebraic
32equations (DAE):
33
34minimize     integral_{t=0}^{10} x0^2 + x1^2 + u^2  dt
35x0,x1,z,u
36
37subject to   dot(x0) == z*x0-x1+u     \
38             dot(x1) == x0             }  for 0 <= t <= 10
39                   0 == x1^2 + z - 1  /
40             x0(t=0) == 0
41             x1(t=0) == 1
42             x0(t=10) == 0
43             x1(t=10) == 0
44             -0.75 <= u <= 1  for 0 <= t <= 10
45
46Note that other methods such as direct collocation or direct multiple shooting
47are usually preferably to the direct single shooting method in practise.
48
49Joel Andersson, 2012-2015
50"""
51
52# Declare variables
53x = SX.sym("x",2) # Differential states
54z = SX.sym("z")   # Algebraic variable
55u = SX.sym("u")   # Control
56
57# Differential equation
58f_x = vertcat(z*x[0]-x[1]+u, x[0])
59
60# Algebraic equation
61f_z = x[1]**2 + z - 1
62
63# Lagrange cost term (quadrature)
64f_q = x[0]**2 + x[1]**2 + u**2
65
66# Create an integrator
67dae = {'x':x, 'z':z, 'p':u, 'ode':f_x, 'alg':f_z, 'quad':f_q}
68opts = {"tf":0.5} # interval length
69I = integrator('I', "idas", dae, opts)
70
71# All controls
72U = MX.sym("U", 20)
73
74# Construct graph of integrator calls
75X  = [0,1]
76J = 0
77for k in range(20):
78  Ik = I(x0=X, p=U[k])
79  X = Ik['xf']
80  J += Ik['qf']   # Sum up quadratures
81
82# Allocate an NLP solver
83nlp = {'x':U, 'f':J, 'g':X}
84opts = {"ipopt.linear_solver":"ma27"}
85solver = nlpsol("solver", "ipopt", nlp, opts)
86
87# Pass bounds, initial guess and solve NLP
88sol = solver(lbx = -0.75, # Lower variable bound
89             ubx =  1.0,  # Upper variable bound
90             lbg =  0.0,  # Lower constraint bound
91             ubg =  0.0,  # Upper constraint bound
92             x0  =  0.0) # Initial guess
93print(sol)
94