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