Remove useless show_pixmap() calls + fix translation
[scilab.git] / scilab / modules / cacsd / demos / pendulum / pendule.dem
1 //
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA
4 //
5 // This file is distributed under the same license as the Scilab package.
6 //
7
8 function demo_pendule()
9
10     mode(-1)
11
12
13     exec("SCI/modules/cacsd/demos/pendulum/simulation.sci", -1);
14     exec("SCI/modules/cacsd/demos/pendulum/graphics.sci", -1);
15
16     //disable displayed lines control
17     display_props=lines(); lines(0)
18     my_handle             = scf(100001);
19     clf(my_handle,"reset");
20
21     show_window();
22     wdim=xget('wdim')
23     //mode(1)
24
25     dpnd()
26
27     //
28     // equations 
29     //----------
30     //state =[x x' theta theta']  
31     //
32     mb=0.1;mc=1;l=0.3;m=4*mc+mb; //constants
33     //
34     messagebox(['Open loop simulation'
35     ' '
36     ' y0 = [0;0;0.1;0]; //initial state [x x'' theta theta'']'
37     ' t = 0.03*(1:180); //observation dates'
38     ' y = ode(y0,0,0.03*(1:180),ivpd); //differential equation integration'
39     ' //Display'
40     ' P = initialize_display(y0(1),y0(3))'
41     ' for k=1:size(y,2), set_pendulum(P,y(1,k),y(3,k));end'],"modal");
42     //
43     y0=[0;0;0.1;0];
44     y=ode(y0,0,0.03*(1:180),ivpd);
45
46     P=initialize_display(y0(1),y0(3));
47     for k=1:size(y,2), set_pendulum(P,y(1,k),y(3,k));end
48     //
49     messagebox(['Linearization'
50     ' '
51     ' x0=[0;0;0;0];u0=0;'
52     ' [f,g,h,j]=lin(pendu,x0,u0);'
53     ' pe=syslin(''c'',f,g,h,j);'
54     ' // display of the linear system'
55     ' ssprint(pe)'],"modal")
56
57     //
58     mode(1)
59     x0=[0;0;0;0];u0=0;
60     [f,g,h,j]=lin(pendu,x0,u0);
61     pe=syslin('c',f,g,h,j);
62     ssprint(pe)
63     mode(-1)
64
65     //
66     messagebox(['Checking the result'
67     ' //comparison with linear model computed by hand';
68     ''
69     ' f1=[0 1        0                0'
70     '     0 0    -3*mb*9.81/m         0'
71     '     0 0        0                1'
72     '     0 0  6*(mb+mc)*9.81/(m*l)   0];'
73     ' '
74     ' g1=[0 ; 4/m ; 0 ; -6/(m*l)];'
75     ' '
76     ' h1=[1 0 0 0'
77     '     0 0 1 0];'
78     ' '
79     ' err=norm(f-f1,1)+norm(g-g1,1)+norm(h-h1,1)+norm(j,1)'],"modal")
80
81
82     //
83     mode(1)
84     f1=[0 1        0             0
85     0 0    -3*mb*9.81/m         0
86     0 0        0             1
87     0 0  6*(mb+mc)*9.81/(m*l)   0];
88     g1=[0 ; 4/m ; 0 ; -6/(m*l)];
89     h1=[1 0 0 0
90     0 0 1 0];
91     err=norm(f-f1,1)+norm(g-g1,1)+norm(h-h1,1)+norm(j,1)
92     mode(-1)
93
94     messagebox(['Linear system properties analysis'
95     ' spec(f) //stability (unstable system !)'
96     ' n=contr(f,g) //controlability'
97     ' '
98     ' //observability'
99     ' m1=contr(f'',h(1,:)'') '
100     ' [m2,t]=contr(f'',h(2,:)'')'],"modal");
101     //---------
102     mode(1)
103     spec(f) //stability (unstable system !)
104     n=contr(f,g) //controlability
105
106     //observability
107     m1=contr(f',h(1,:)') 
108     [m2,t]=contr(f',h(2,:)')
109     mode(-1)
110     //
111     messagebox(['Synthesis of a stabilizing controller using poles placement'
112     ' '
113     '// only x and theta are observed  : contruction of an observer'
114     '// to estimate the state : z''=(f-k*h)*z+k*y+g*u'
115     ' to=0.1; ' 
116     ' k=ppol(f'',h'',-ones(4,1)/to)''  //observer gain'
117     '// compute the compensator gain'
118     ' kr=ppol(f,g,-ones(4,1)/to)'],"modal");
119
120     //-------------------------------------------------
121     //
122     //pole placement technique
123     //only x and theta are observed  : contruction of an observer
124     //to estimate the state : z'=(f-k*h)*z+k*y+g*u
125     //
126     to=0.1;  //
127     k=ppol(f',h',-ones(4,1)/to)'  //observer gain
128     //
129     //verification
130     //
131     // norm( poly(f-k*h,'z')-poly(-ones(4,1)/to,'z'))
132     //
133     kr=ppol(f,g,-ones(4,1)/to)  //compensator gain
134
135     //
136     messagebox(['Build full linear system  pendulum-observer-compensator'
137     ' '
138     ' ft=[f-g*kr            -g*kr'
139     '     0*f               f-k*h]'
140     ' '
141     ' gt=[g;0*g];'
142     ' ht=[h,0*h];'
143     ''
144     ' pr=syslin(''c'',ft,gt,ht);'
145     ''
146     '//Check the closed loop dynamic'
147     ' spec(pr.A)'
148     ' '
149     '//transfer matrix representation'
150     ' hr=clean(ss2tf(pr),1.d-10)'
151     ' '
152     '//frequency analysis'
153     ' black(pr,0.01,100,[''position'',''theta''])'
154     ' g_margin(pr(1,1))'],"modal");
155
156     //---------------------------------------------
157     //
158     //state: [x x-z]
159     //
160     mode(1)
161     ft=[f-g*kr            -g*kr
162     0*f               f-k*h];
163     gt=[g;0*g];
164     ht=[h,0*h];
165     pr=syslin('c',ft,gt,ht);
166
167     // closed loop dynamics:
168     spec(pr(2))
169     //transfer matrix representation
170     hr=clean(ss2tf(pr),1.d-10)
171
172
173     //frequency analysis
174     clf()
175     black(pr,0.01,100,['position','theta'])
176     g_margin(pr(1,1))
177     mode(-1)
178
179
180     //
181     messagebox(['Sampled system'
182     ' '
183     ' t=to/5; //sampling period'
184     ' prd=dscr(pr,t); //discrete model'
185     ' spec(prd.A) //poles of the discrete model'],"modal");
186
187     //---------------
188     //
189     mode(1)
190     t=to/5;
191     prd=dscr(pr,t);
192     spec(prd(2))
193     mode(-1)
194     //
195     messagebox(['Impulse response'
196     ' '
197     ' x0=[0;0;0;0;0;0;0;0]; //initial state'
198     ' u(1,180)=0;u(1,1)=1; //impulse'
199     ' y=flts(u,prd,x0);    //siscrete system simulation'
200     ' draw1()'],"modal");
201
202     //-----------------
203     //
204     mode(1)
205     x0=[0;0;0;0;0;0;0;0];
206     u(1,180)=0;u(1,1)=1;
207     y=flts(u,prd,x0);
208
209     draw1()
210     mode(-1)
211     //
212     messagebox(['Compensation of the non linear system with'
213     'linear regulator'
214     ''
215     ' t0=0;t1=t*(1:125);'
216     ' x0=[0 0 0.4 0   0 0 0 0]'';  // initial state'
217     ' yd=ode(x0,t0,t1,regu); //integration of differential equation'
218     ' draw2()'],"modal");
219     ;
220     //--------------------------------------
221     //
222     //simulation
223     //
224     mode(1)
225     t0=0;t1=t*(1:125);
226     x0=[0 0 0.4 0   0 0 0 0]';   //
227     yd=ode(x0,t0,t1,regu);
228     draw2()
229     mode(-1) 
230     lines(display_props(2))
231     messagebox('The end',"modal");
232 endfunction
233
234 demo_pendule();
235 clear demo_pendule;