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