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