/home/casadibot/slaves/linuxbot/full/source/docs/tutorials/python/src/integration/temp.py
January 17, 2015
CasADi tutorial
0
1
2
3
4
5
6
11
12
13
14
#
#
#
#
#
#
#
1
There are two basic mechanisms two retrieve the whole trajectory of states:
- method A, using reset + integrate
- method B, using Simulator
We demonstrate first method A:
from numpy i m p o r t *
i m p o r t numpy
from c a s a d i i m p o r t *
from p y l a b i m p o r t *
47
48
49
50
51
52
t s =numpy . l i n s p a c e (0, t e n d ,100)
x0 = 0; y0 = 1
i n t e g r a t o r . s e t I n p u t ( [ x0 , y0 ] ," x0 ")
i n t e g r a t o r . s e t I n p u t (0," p ")
i n t e g r a t o r . e v a l u a t e ()
i n t e g r a t o r . r e s e t ()
Define a convenience function to acces x(t)
57
def out( t ) :
Error:
58
integrator.integrate(t)
59
r e t u r n i n t e g r a t o r . o u t p u t (). t o A r r a y ()
/usr/lib/pymodules/python2.7/matplotlib/__init__.py:758: UserWarning: Found matplotlib
configuration in ~/.matplotlib/. To conform with the XDG base directory standard,
60
_get_xdg_config_dir())
61
s o l = array( [ o u t ( t ) f o r t i n t s ] ). s q u e e z e ()
Plot the trajectory
ODE integration
Let’s construct a simple Van der Pol oscillator.
18
19
20
21
25
26
u = SX. sym (" u ")
x = SX. sym (" x ")
y = SX. sym (" y ")
f
= S X F u n c t i o n ( [ v e r t c a t ( [ x , y ] ), u ] , [ v e r t c a t ( [ (1- y * y )* x - y + u , x ] ) ] )
62
63
64
65
66
67
f i g u r e ()
p l o t ( s o l [ : ,0 ] , s o l [ : ,1 ] )
t i t l e ( ’ Van d e r P o l p h a s e s p a c e ’)
x l a b e l ( ’ x ’)
y l a b e l ( ’ y ’)
show ()
Manipulate the function to adhere to the integrator’s input/output signature f(time;states;parameters)
2.0
fmod = S X F u n c t i o n ( d a e I n ( x = f . i n p u t E x p r (0), p= f . i n p u t E x p r (1)), d a e O u t ( ode = f .
o u t p u t E x p r (0)))
fmod . s e t O p t i o n (" name ","ODE r i g h t hand s i d e ")
1.5
1.0
Create the Integrator
28
i n t e g r a t o r = I n t e g r a t o r (" c v o d e s ", fmod )
0.5
The whole series of sundials options are available for the user
0.0
i n t e g r a t o r . s e t O p t i o n (" f s e n s _ e r r _ c o n ", T r u e )
i n t e g r a t o r . s e t O p t i o n (" q u a d _ e r r _ c o n ", T r u e )
i n t e g r a t o r . s e t O p t i o n (" a b s t o l ",1 e -6)
i n t e g r a t o r . s e t O p t i o n (" r e l t o l ",1 e -6)
t e n d = 10
i n t e g r a t o r . s e t O p t i o n (" t 0 ",0)
i n t e g r a t o r . s e t O p t i o n (" t f ", t e n d )
i n t e g r a t o r . i n i t ()
p r i n t isinstance( i n t e g r a t o r , F u n c t i o n )
y
30
31
32
33
34
35
36
37
38
0.5
1.0
1.5
2.0
True
39
2.5
p r i n t "%d −> %d " % ( i n t e g r a t o r . g e t N u m I n p u t s (), i n t e g r a t o r . g e t N u m O u t p u t s ())
6 -> 6
Setup the Integrator to integrate from 0 to t=tend, starting at [x0,y0]
The output of Integrator is the state at the end of integration.
Van der Pol phase space
3
2
We demonstrate method B:
70
71
72
sim = S i m u l a t o r ( i n t e g r a t o r , t s )
s i m . i n i t ()
s i m . s e t I n p u t ( [ x0 , y0 ] ," x0 ")
1
0
x
1
2
3
/home/casadibot/slaves/linuxbot/full/source/docs/tutorials/python/src/integration/temp.py
73
74
75
76
s i m . s e t I n p u t (0," p ")
s i m . e v a l u a t e ()
78
print
99
100
101
102
103
104
105
106
107
108
109
s o l 2 = s i m . o u t p u t (). t o A r r a y ().T
sol and sol2 are exactly the same
l i n a l g . norm ( s o l - s o l 2 )
0.0
Sensitivity for initial conditions
2
#
d i n t e g r a t o r = i n t e g r a t o r . d e r i v a t i v e (1,0)
d i n t e g r a t o r . s e t I n p u t ( [ x0 , y0 ] ," d e r _ x 0 ")
d i n t e g r a t o r . s e t I n p u t ( [ 1,0 ] ," fwd0_x0 ")
d i n t e g r a t o r . e v a l u a t e ()
A = d i n t e g r a t o r . g e t O u t p u t (" f w d 0 _ x f ") [ 0 ]
p l o t ( dx0 ,A* dx0 )
l e g e n d (( ’ T r u e s e n s i t i v i t y ’, ’ L i n e a r i s e d s e n s i t i v i t y ’))
p l o t (0,0, ’ o ’)
show ()
0.5
0.0
o u t = array( [ o u t ( dx ) f o r dx i n dx0 ] ). s q u e e z e ()
0.5
dx(tend)
d x t e n d = o u t [ : ,0 ] - s o l [ -1,0 ]
f i g u r e ()
p l o t ( dx0 , d x t e n d )
g r i d ()
t i t l e ( ’ I n i t i a l p e r t u r b a t i o n map ’)
x l a b e l ( ’ dx ( 0 ) ’)
y l a b e l ( ’ dx ( t e n d ) ’)
show ()
2.0
2.5
3.0
2.0
0.0
1.5
1.0
0.5
0.0
dx(0)
0.5
1.0
1.5
2.0
The interpetation is that a small initial circular patch of phase space evolves into ellipsoid patches at later stages.
0.5
1.0
1.5
2.0
2.5
3.0
2.0
1.0
1.5
Initial perturbation map
0.5
Initial perturbation map
True sensitivity
Linearised sensitivity
1.0
d e f o u t ( dx0 ) :
i n t e g r a t o r . s e t I n p u t ( [ x0 + dx0 , y0 ] ," x0 ")
i n t e g r a t o r . e v a l u a t e ()
r e t u r n i n t e g r a t o r . o u t p u t (). t o A r r a y ()
dx0 =numpy . l i n s p a c e (-2,2,100)
dx(tend)
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
January 17, 2015
1.5
1.0
0.5
0.0
dx(0)
0.5
1.0
1.5
2.0
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def out( t ) :
d i n t e g r a t o r . s e t I n p u t ( [ 1,0 ] ," fwd0_x0 ")
d i n t e g r a t o r . e v a l u a t e ()
A= d i n t e g r a t o r . o u t p u t (" f w d 0 _ x f "). t o A r r a y ()
d i n t e g r a t o r . s e t I n p u t ( [ 0,1 ] ," fwd0_x0 ")
d i n t e g r a t o r . e v a l u a t e ()
B= d i n t e g r a t o r . o u t p u t (" f w d 0 _ x f "). t o A r r a y ()
r e t u r n array( [ A,B ] ). s q u e e z e ().T
c i r c l e = array( [ [ s i n ( x ), c o s ( x ) ]
f o r x i n numpy . l i n s p a c e (0,2* p i ,100) ] ).T
f i g u r e ()
p l o t ( s o l [ : ,0 ] , s o l [ : ,1 ] )
g r i d ()
f o r i i n range(10) :
J = o u t ( t s [ 10* i ] )
e= 0.1* numpy . d o t ( J , c i r c l e ).T+ s o l [ 10* i , : ]
/home/casadibot/slaves/linuxbot/full/source/docs/tutorials/python/src/integration/temp.py
138
139
140
January 17, 2015
p l o t ( e [ : ,0 ] , e [ : ,1 ] , c o l o r = ’ r e d ’)
3
[ 1.39637624 ] ]
show ()
3
2
1
0
1
Let’s plot the results
2
3
3
2
1
0
1
2
3
The figure reveals that perturbations perpendicular to the phase space trajectory shrink.
Symbolic intergator results
164
165
166
167
168
169
170
171
172
173
174
uv =numpy . l i n s p a c e (-1,1,100)
o u t = array( [ o u t ( i ) f o r i i n uv ] ). s q u e e z e ()
f i g u r e ()
p l o t ( uv , o u t )
g r i d ()
t i t l e ( ’ Dependence o f f i n a l s t a t e on i n p u t u ’)
x l a b e l ( ’ u ’)
y l a b e l ( ’ s t a t e ’)
l e g e n d (( ’ x ’, ’ y ’))
show ()
Since Integrator is just another Function, the usual CasADi rules for symbolic evaluation are active.
We create an MX ’w’ that contains the result of a time integration with: - a fixed integration start time , t=0s
- a fixed integration end time, t=10s - a fixed initial condition (1,0) - a free symbolic input, held constant during
integration interval
148
149
u=MX. sym (" u ")
w, = i n t e g r a t o r O u t ( i n t e g r a t o r . c a l l ( i n t e g r a t o r I n ( x0 =MX( [ 1,0 ] ), p=u ))," x f ")
152
153
154
155
156
157
158
159
160
f = MXFunction( [ u ] , [w] )
f . i n i t ()
2
print
o u t (0)
x
y
1
We construct an MXfunction and a python help function ’out’
0
state
d e f o u t (u) :
f . s e t I n p u t (u)
f . e v a l u a t e ()
r e t u r n f . o u t p u t (). t o A r r a y ()
Dependence of final state on input u
1
2
3
[ [ -2.54395395 ]
[ -0.43932676 ] ]
161
print
o u t (1)
[ [ -0.25397819 ]
4
1.0
0.5
0.0
u
0.5
1.0
© Copyright 2024 ExpyDoc
