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