integrators/idas.py

This example looks at the use of casadi::Integrator();

View output (PDF) | source (python)

See also
CASADI_EXPORT Function integrator(const std::string &name, const std::string &solver, const SXDict &dae, const Dict &opts=Dict())
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 from casadi import *
32 from numpy import *
33 from pylab import *
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 L = SX.sym("L")
44 g = SX.sym("g")
45 
46 
47 
48 x=SX.sym("x")
49 y=SX.sym("y")
50 u=SX.sym("u")
51 v=SX.sym("v")
52 
53 
54 
55 lambd=SX.sym("lambda")
56 
57 
58 
59 x_all = vertcat(x,u,y,v)
60 z_all = lambd
61 p_all = vertcat(L,g)
62 
63 
64 
65 P_ = [5,10]
66 X_ = [3,-1.0/3,4,1.0/4]
67 XDOT_ = [-1.0/3,1147.0/240,1.0/4,-653.0/180]
68 Z_ = [1147.0/720]
69 
70 
71 
72 ode = vertcat(u,lambd*x,v,lambd*y+g)
73 alg = x**2+y**2-L**2
74 dae = {'x':x_all, 'z':z_all, 'p':p_all, 'ode':ode, 'alg':alg}
75 f = Function('f', [x_all, z_all, p_all], [ode, alg], ['x', 'z', 'p'], ['ode', 'alg'])
76 
77 
78 
79 res = f(p=P_, x=X_, z=Z_)
80 print(res['ode'])
81 print(res['alg'])
82 
83 
84 
85 
86 j = jacobian(alg,lambd)
87 print(j)
88 
89 
90 
91 
92 
93 
94 
95 
96 I = integrator('I', 'idas', dae, {'calc_ic':False, 'init_xdot':XDOT_})
97 
98 
99 
100 
101 try:
102  I(p=P_, x0=X_, z0=Z_)
103 except Exception as e:
104  print(e)
105 
106 
107 
108 ode = vertcat(u,lambd*x)
109 alg = vertcat(x**2+y**2-L**2, u*x+v*y,u**2-g*y+v**2+L**2*lambd)
110 x_all = vertcat(x,u)
111 z_all = vertcat(y,v,lambd)
112 dae = {'x':x_all, 'z':z_all, 'p':p_all, 'ode':ode, 'alg':alg}
113 f = Function('f', [x_all, z_all, p_all], [ode, alg], ['x', 'z', 'p'], ['ode', 'alg'])
114 
115 
116 
117 P_ = [5,10]
118 X_ = [3,-1.0/3]
119 XDOT_ = [-1.0/3,1147.0/240]
120 Z_ = [4,1.0/4,1147.0/720]
121 
122 
123 
124 res = f(p=P_, x=X_, z=Z_)
125 print(res['ode'])
126 print(res['alg'])
127 
128 
129 
130 J = f.factory('J', f.name_in(), ['jac:alg:z'])
131 res = J(p=P_, x=X_, z=Z_)
132 print(array(res["jac_alg_z"]))
133 
134 
135 
136 
137 
138 I = integrator('I', 'idas', dae, {'t0':0, 'tf':1, 'init_xdot':XDOT_})
139 res = I(p=P_, x0=X_, z0=Z_)
140 print(res['xf'])
141 
142 
143 
144 
145 
146 
147 P_ = [5,10]
148 X_ = [5,0]
149 
150 
151 
152 try:
153  I(p=P_, x0=X_, z0=Z_)
154 except Exception as e:
155  print(e)
156 
157 
158