bug 7514 fixed - Some demo files were broken due to changes in x_mdialog function.
[scilab.git] / scilab / modules / cacsd / demos / pid.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 mode(-1);
9 lines(0);
10 //display the diagram
11 x=[5,10,20,40,50,70,80,90];xmin=-10;xmax=100;
12 y=[22,28,30,32];ymin=12;ymax=40;
13
14 xx=[xmin,xmin,x([1 2 2 7 4 6 3 4 5 6 3 3 5 5]);xmax,xmax,x([3,2,7,7,5,8,3,4,5,6,4,4,6,6])];
15 yy=[ymin,ymax,y([3,1,1,1,3,3,2,2,2,2,2,4,2,4]);ymin,ymax,y([3,3,1,3,3,3,4,4,4,4,2,4,2,4])];
16
17 xset("window",0);clf();show_window();
18 plot2d(xx,yy,ones(1,16),'022');
19 xstring(28,30,'K');xstring(56,30,'Plant');xstring(12,28.80,'-');
20 xtitle('PLANT   and   CONTROLLER')
21 mode(2);
22
23 s=poly(0,'s');z=poly(0,'z');
24 messagebox(['Example of PID Design '
25             'file: '+SCI+'/modules/cacsd/demos/pid.dem'],"modal");
26
27 n=x_choose(['Continuous time';'Discrete time'],'Select time domain');
28 select n
29  case 0
30   warning('Demo stops!');return;
31  case 1
32   mode(1)
33   dom='c';
34   s=poly(0,'s');
35   str='[(s-1)/(s^2+5*s+1)]';
36   rep=x_dialog('Nominal plant?',str);
37   if rep==[] then return,end
38   Plant=evstr(rep);
39   Plant=syslin('c',Plant);
40   mode(-1)
41 case 2
42   mode(1)
43   dom='d'
44   z=poly(0,'z');
45   str='(z+1)/(z^2-5*z+2)'
46   rep=x_dialog('Nominal plant?',str);
47   if rep==[] then return,end
48   Plant=evstr(rep)
49   Plant=syslin('d',Plant);
50   mode(-1)
51 end
52 //Nominal Plant
53 P22=tf2ss(Plant);    //...in state-space form
54 [ny,nu,nx]=size(P22);
55 defv=['-1.2';'1';'0.1'];
56 if dom=='d' then defv=['-10';'1';'0.1'];end
57 while %t
58   mode(1)
59   if dom=='c' then
60     _title='Enter your PID controller K(s)=Kp*(1+T0/s+T1*s)';
61   end
62   if dom=='d' then
63     _title='Enter your PID controller K(z)=Kp*(1+T0/z+T1*z)';
64   end
65   defv=x_mdialog(_title,['Kp=';'T0=';'T1='],defv);
66   if defv==[] then warning('Demo stops!');return;end
67   Kp=evstr(defv(1));T0=evstr(defv(2));T1=evstr(defv(3));
68   if dom=='c' then
69     Kpid=tf2ss(Kp*(1+T0/s+T1*s));
70   end
71   if dom=='d' then
72     Kpid=tf2ss(Kp*(1+T0/z+T1*z));
73   end
74   W=[1, -P22;
75       Kpid,1];Winv=inv(W);
76
77   disp(spec(Winv(2)),'closed loop eigenvalues');//Check internal stability
78   if max(real(spec(Winv(2)))) > 0 then 
79     messagebox('You loose: closed-loop is UNSTABLE!!!',"modal");
80   else
81     messagebox('Congratulations: closed-loop is STABLE !!!',"modal");
82     break;
83   end
84   mode(-1)
85 end
86 mode(1)
87 [Spid,Rpid,Tpid]=sensi(P22,Kpid);  //Sensitivity functions
88 Tpid(5)=clean(Tpid(5));
89
90 disp(clean(ss2tf(Spid)),'Sensitivity function');
91 disp(clean(ss2tf(Tpid)),'Complementary sensitivity function');
92
93 resp=['Frequency response';'Time response'];
94 while %t do
95   n=x_choose(resp,'Select response(s)');
96   if degree(Tpid(5))>0 then
97     warning('Improper transfer function! D(s) set to D(0)')
98     Tpid(5)=coeff(Tpid(5),0);
99   end
100   Tpid(5)=coeff(Tpid(5));
101   select n
102   case 0
103     break
104   case 1
105     mode(1)
106     clf(1);xset("window",1);show_window();bode(Tpid);
107     mode(-1)
108   case 2
109     if Plant(4)=='c' then
110       mode(1)
111       defv=['0.1','50'];
112       _title='Enter Sampling period and Tmax';
113       rep=x_mdialog(_title,['Sampling period?';'Tmax?'],defv);
114       if rep==[] then break,end
115       dttmax=evstr(rep);
116       dt=evstr(dttmax(1));tmax=evstr(dttmax(2));
117       t=0:dt/5:tmax;
118       n1=x_choose(['Step response?';'Impulse response?'],'Simulation:');
119       if n1==0 then
120         warning('Demo stops!');return;
121       end
122       if n1==1 then 
123         clf(1);xset("window",1);show_window();
124         plot2d([t',t'],[(csim('step',t,Tpid))',ones(t')])
125       end
126       if n1==2 then
127         clf(1);xset("window",1);show_window();
128         plot2d([t',t'],[(csim('impul',t,Tpid))',0*t'])
129       end
130       mode(-1)
131     elseif Plant(4)=='d' then
132       mode(1)
133       defv=['100'];
134       _title='Tmax?'
135       rep=x_mdialog(_title,['Tmax='],defv);
136       if rep==[] then break,end
137       Tmax=evstr(rep);
138       mode(-1)
139       while %t do
140         n=x_choose(['Step response?';'Impulse response?'],'Simulation:');
141         select n
142         case 0 then
143           break
144         case 1 then
145           mode(1)
146           u=ones(1,Tmax);u(1)=0;
147           clf(1);xset("window",1);show_window();
148           plot2d([(1:Tmax)',(1:Tmax)'],[(dsimul(Tpid,u))',(ones(1:Tmax)')])
149           mode(-1)
150         case 2 then
151           mode(1)
152           u=zeros(1,Tmax);u(1)=1;
153           clf(1);xset("window",1);show_window();
154           plot2d((1:Tmax)',(dsimul(Tpid,u))')
155           mode(-1)
156         end
157       end
158     end
159   end
160 end