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#
24# -*- coding: utf-8 -*-
25
26from __future__ import print_function
27from casadi import *
28from copy import deepcopy
29print('Testing sensitivity analysis in CasADi')
30
31# All ODE and DAE integrators to be tested
32DAE_integrators = ['idas','collocation']
33ODE_integrators = ['cvodes','rk'] + DAE_integrators
34
35for Integrators in (ODE_integrators,DAE_integrators):
36  if Integrators==ODE_integrators: # rocket example
37    print('******')
38    print('Testing ODE example')
39
40    # Time
41    t = SX.sym('t')
42
43    # Parameter
44    u = SX.sym('u')
45
46    # Differential states
47    s = SX.sym('s'); v = SX.sym('v'); m = SX.sym('m')
48    x = vertcat(s,v,m)
49
50    # Constants
51    alpha = 0.05 # friction
52    beta = 0.1   # fuel consumption rate
53
54    # Differential equation
55    ode = vertcat(
56      v,
57      (u-alpha*v*v)/m,
58      -beta*u*u)
59
60    # Quadrature
61    quad = v**3 + ((3-sin(t)) - u)**2
62
63    # DAE callback function
64    dae = {'t':t, 'x':x, 'p':u, 'ode':ode, 'quad':quad}
65
66    # Time length
67    tf = 0.5
68
69    # Initial position
70    x0 = [0.,0.,1.]
71
72    # Parameter
73    u0 = 0.4
74
75  else: # Simple DAE example
76    print('******')
77    print('Testing DAE example')
78
79    # Differential state
80    x = SX.sym('x')
81
82    # Algebraic variable
83    z = SX.sym('z')
84
85    # Parameter
86    u = SX.sym('u')
87
88    # Differential equation
89    ode = -x + 0.5*x*x + u + 0.5*z
90
91    # Algebraic constraint
92    alg = z + exp(z) - 1.0 + x
93
94    # Quadrature
95    quad = x*x + 3.0*u*u
96
97    # DAE callback function
98    dae = {'x':x, 'z':z, 'p':u, 'ode':ode, 'alg':alg, 'quad':quad}
99
100    # End time
101    tf = 5.
102
103    # Initial position
104    x0 = 1.
105
106    # Parameter
107    u0 = 0.4
108
109  # Integrator
110  for MyIntegrator in Integrators:
111    print('========')
112    print('Integrator: ', MyIntegrator)
113    print('========')
114
115    # Integrator options
116    opts = {'tf':tf}
117    if MyIntegrator=='collocation':
118      opts['rootfinder'] = 'kinsol'
119
120    # Integrator
121    I = integrator('I', MyIntegrator, dae, opts)
122
123    # Integrate to get results
124    res = I(x0=x0, p=u0)
125    xf = res['xf']
126    qf = res['qf']
127    print('%50s' % 'Unperturbed solution:', 'xf  = ', xf, ', qf  = ', qf)
128
129    # Perturb solution to get a finite difference approximation
130    h = 0.001
131    res = I(x0=x0, p=u0+h)
132    fd_xf = (res['xf']-xf)/h
133    fd_qf = (res['qf']-qf)/h
134    print('%50s' % 'Finite difference approximation:', 'd(xf)/d(p) = ', fd_xf, ', d(qf)/d(p) = ', fd_qf)
135
136    # Calculate once, forward
137    I_fwd = I.factory('I_fwd', ['x0', 'z0', 'p', 'fwd:p'], ['fwd:xf', 'fwd:qf'])
138    res = I_fwd(x0=x0, p=u0, fwd_p=1)
139    fwd_xf = res['fwd_xf']
140    fwd_qf = res['fwd_qf']
141    print('%50s' % 'Forward sensitivities:', 'd(xf)/d(p) = ', fwd_xf, ', d(qf)/d(p) = ', fwd_qf)
142
143    # Calculate once, adjoint
144    I_adj = I.factory('I_adj', ['x0', 'z0', 'p', 'adj:qf'], ['adj:x0', 'adj:p'])
145    res = I_adj(x0=x0, p=u0, adj_qf=1)
146    adj_x0 = res['adj_x0']
147    adj_p = res['adj_p']
148    print('%50s' % 'Adjoint sensitivities:', 'd(qf)/d(x0) = ', adj_x0, ', d(qf)/d(p) = ', adj_p)
149
150    # Perturb adjoint solution to get a finite difference approximation of the second order sensitivities
151    res = I_adj(x0=x0, p=u0+h, adj_qf=1)
152    fd_adj_x0 = (res['adj_x0']-adj_x0)/h
153    fd_adj_p = (res['adj_p']-adj_p)/h
154    print('%50s' % 'FD of adjoint sensitivities:', 'd2(qf)/d(x0)d(p) = ', fd_adj_x0, ', d2(qf)/d(p)d(p) = ', fd_adj_p)
155
156    # Forward over adjoint to get the second order sensitivities
157    I_foa = I_adj.factory('I_foa', ['x0', 'z0', 'p', 'adj_qf', 'fwd:p'], ['fwd:adj_x0', 'fwd:adj_p'])
158    res = I_foa(x0=x0, p=u0, adj_qf=1, fwd_p=1)
159    fwd_adj_x0 = res['fwd_adj_x0']
160    fwd_adj_p = res['fwd_adj_p']
161    print('%50s' % 'Forward over adjoint sensitivities:', 'd2(qf)/d(x0)d(p) = ', fwd_adj_x0, ', d2(qf)/d(p)d(p) = ', fwd_adj_p)
162
163    # Adjoint over adjoint to get the second order sensitivities
164    I_aoa = I_adj.factory('I_aoa', ['x0', 'z0', 'p', 'adj_qf', 'adj:adj_p'], ['adj:x0', 'adj:p'])
165    res = I_aoa(x0=x0, p=u0, adj_qf=1, adj_adj_p=1)
166    adj_x0 = res['adj_x0']
167    adj_p = res['adj_p']
168    print('%50s' % 'Adjoint over adjoint sensitivities:', 'd2(qf)/d(x0)d(p) = ', adj_x0, ', d2(qf)/d(p)d(p) = ', adj_p)
169