← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to demos/control/lqg.dem

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
demos/control/lqg.dem
 
1
// Copyright INRIA
 
2
mode(-1)
 
3
exec(SCI+'/demos/control/scheme.dem');
 
4
s=poly(0,'s');z=poly(0,'z');
 
5
x_message(['Simple example of SISO LQG Design';
 
6
           'Demo is in file '+SCI+'/demos/control/lqg.dem';
 
7
           'Computes the LQG compensator and plots response'])
 
8
n=x_choose(['Continuous time';'Discrete time'],'Select time domain');
 
9
//  
 
10
select n
 
11
case 0
 
12
  return
 
13
case 1
 
14
  mode(1)
 
15
  s=poly(0,'s');
 
16
  str='(s-1)/(s^2-5*s+1)';
 
17
  rep=x_dialog('Nominal plant?',str)
 
18
  if rep==[] then return,end
 
19
  Plant=evstr(rep);
 
20
  Plant=syslin('c',Plant);
 
21
  mode(-1)
 
22
case 2
 
23
  mode(1)
 
24
  z=poly(0,'z');
 
25
  str='(z+1)/(z^2-5*z+2)'
 
26
  rep=x_dialog('Nominal plant?',str)
 
27
  if rep==[] then return,end
 
28
  Plant=evstr(rep);
 
29
  Plant=syslin('d',Plant);
 
30
  mode(-1)
 
31
end   
 
32
mode(1)
 
33
//Nominal Plant
 
34
P22=tf2ss(Plant);    //...in state-space form
 
35
[ny,nu,nx]=size(P22);
 
36
x_message('Now enter weighting matrices');
 
37
rep=x_matrix('x-weighting matrix',eye(nx,nx))
 
38
if rep==[] then return,end
 
39
Qx=evstr(rep);
 
40
rep=x_matrix('u-weighting matrix',eye(nu,nu));
 
41
if rep==[] then return,end
 
42
Qu=evstr(rep);
 
43
bigQ=sysdiag(Qx,Qu);
 
44
rep=x_matrix('x-noise covariance matrix',eye(nx,nx))
 
45
if rep==[] then return,end
 
46
Rx=evstr(rep);
 
47
rep=x_matrix('y-noise covariance matrix',eye(ny,ny))
 
48
if rep==[] then return,end
 
49
Ry=evstr(rep);
 
50
bigR=sysdiag(Rx,Ry);
 
51
[Plqg,r]=lqg2stan(P22,bigQ,bigR);     //LQG pb as a standard problem
 
52
Klqg=lqg(Plqg,r);                     //LQG compensator
 
53
 
 
54
disp(spec(h_cl(Plqg,r,Klqg)),'closed loop eigenvalues:');    //Check internal stability
 
55
[Slqg,Rlqg,Tlqg]=sensi(P22,Klqg);  //Sensitivity functions
 
56
 
 
57
disp(clean(ss2tf(Slqg)),'Sensitivity function');
 
58
disp(clean(ss2tf(Tlqg)),'Complementary sensitivity function');
 
59
mode(-1)
 
60
x_message('Closed-loop response');
 
61
 
 
62
resp=['Frequency response';'Time response'];
 
63
while %t do
 
64
  n=x_choose(resp,'Select response(s)');
 
65
  select n
 
66
  case 0
 
67
    disp("LQG demo stops!");break;
 
68
  case 1 then
 
69
    mode(1)
 
70
    xbasc(0);xset("window",0);xselect();bode(Tlqg);xend()
 
71
    mode(-1)
 
72
  case 2
 
73
    if Plant(4)=='c' then
 
74
      mode(1)
 
75
      defv=['0.1','20'];
 
76
      title='Enter Sampling period and Tmax';
 
77
      rep=x_mdialog(title,['Sampling period?';'Tmax?'],defv)
 
78
      if rep==[] then break,end
 
79
      dttmax=evstr(rep)
 
80
      dt=evstr(dttmax(1));tmax=evstr(dttmax(2));
 
81
      t=0:dt/5:tmax;
 
82
      n1=x_choose(['Step response?';'Impulse response?'],'Simulation:');
 
83
      select n1
 
84
      case 1 then
 
85
        xbasc(1);xset("window",1);xselect();
 
86
        plot2d([t',t'],[(csim('step',t,Tlqg))',ones(t')]);xend()
 
87
      case 2 then
 
88
        xbasc(1);xset("window",1);xselect();
 
89
        plot2d([t',t'],[(csim('impul',t,Tlqg))',0*t']);xend()
 
90
      end
 
91
      mode(-1)
 
92
    elseif Plant(4)=='d' then
 
93
      mode(1)
 
94
      defv=['30'];
 
95
      title='Tmax?';
 
96
      rep=x_mdialog(title,['Tmax='],defv)
 
97
      if rep==[] then break,end
 
98
      Tmax=evstr(rep);
 
99
      mode(-1)
 
100
      n2=x_choose(['Step response?';'Impulse response?'],'Simulation:');
 
101
      select n2
 
102
      case 0 then
 
103
        break;
 
104
      case 1 then
 
105
        mode(1)
 
106
        u=ones(1,Tmax);u(1)=0;
 
107
        xbasc(1);xset("window",1);xselect();
 
108
        plot2d([(1:Tmax)',(1:Tmax)'],[(dsimul(Tlqg,u))',(ones(1:Tmax)')])
 
109
        mode(-1)
 
110
        xend()  
 
111
      case 2 then
 
112
        mode(1)
 
113
        u=zeros(1,Tmax);u(1)=1;
 
114
        xbasc(1);xset("window",1);xselect();
 
115
        plot2d((1:Tmax)',(dsimul(Tlqg,u))')
 
116
        mode(-1)
 
117
        xend()  
 
118
      end
 
119
    end
 
120
  end
 
121
end
 
122
clear s z n str rep Plant P22 ny nu nx Qx Qu bigQ Rx Ry bigR  Plqg r 
 
123
clear Klqg Slqg Rlqg Tlqg resp title dttmax dt t n1 defv Tmax n2 u