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