ep10_1.m

离散控制系统设计的MATLAB 代码

%%%%%%%%%%%% Exploratory problem 10.1 %%%%%%%%%%%
%   Discrete-Time Control Problems using        %
%       MATLAB and the Control System Toolbox   %
%   by J.H. Chow, D.K. Frederick, & N.W. Chbat  %
%         Brooks/Cole Publishing Company        %
%                September 2002                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ---- Interactive analysis for state-space models ----
%
clear,clf
activity = -1;      % force activity menu to appear
while (activity < 1) | (activity > 7),
  disp(' 1 <= Enter state-space models')
  disp(' 2 <= Check controllability and observability')
  disp(' 3 <= Design full-state feedback gain using pole-placement')
  disp(' 4 <= Design observer gain using pole-placement')
  disp(' 5 <= Form observer controller and see its frequency response')
  disp(' 6 <= Initial-condition response')
  disp(' 7 <= Step response')
  activity = input('Select action, or enter 0 to quit ==> ');

if activity == 0 
  disp('Completed EP10.1 on state-space model design'), return
end

if activity == 1,	% Enter state-space models
model_entry = -1;   % force model entry menu to appear
while (model_entry < 1) | (model_entry > 6),
  disp(' 1 <= Enter A,B,C,D matrices and sampling period')
  disp(' 2 <= Activate keyboard to load a data file')
  model_entry = input...
    ('Select model entry or enter 0 for top-level menu ==> '); 

if (model_entry == 0 | isempty(model_entry)),
  activity = -1;    % return to activity menu
  break             % terminate "while model_entry..." loop
end   % of " if model_entry == 0"

%+++++++++++++++++++ enter model in SS form ++++++++++++++++++++
  if model_entry == 1,
    disp(' ')
  disp('   **** Interactive entry of transfer function in SS form ****')
  disp(' ')
% remove any stale zeros, poles, & gain from workspace
  disp('Enter coefficients of A,B,C,D matrices and sampling period Ts.')
  disp('Have initialized with values from Example 10.6')
  %-----------------------------------------------------
  A_last = [0.81   -0.23   -0.045       % A matrix in Example 10.6
            0.09    0.98   -0.0023	
            0.005   0.10    1     ];
  B_last = [0.09 0.0047 0.00016]';      % B matrix
  C_last = [1 3.5 3];                   % C matrix
  D_last = 0;                           % D matrix
  Ts_last = 0.1;
  cmd = 'n';	% 'n' ==> model is not OK
  while cmd == 'n',
    disp(' ')
    % allow user to enter new values for matrices
    A = input('Enter A matrix ==> ');
    if isempty(A), A = A_last; end      % CR ==> use old value
    B = input('Enter B matrix ==> ');
    if isempty(B), B = B_last; end      % CR ==> use old value
    C = input('Enter C matrix ==> ');
    if isempty(C), C = C_last; end      % CR ==> use old value
    D = input('Enter D matrix ==> ');
    if isempty(D), D = D_last; end      % CR ==> use old value
    Ts = input('Enter sampling period ==> ');
    if isempty(Ts), Ts = Ts_last; end   % CR ==> use old value
    disp(' ')
    disp('A,B,C,D matrices of SS form are:')
    A, B, C, D, Ts
    cmd = input('Model OK? (y|n)','s');
    A_last = A;    % save previous values in case want to reuse
    B_last = B;
    C_last = C;
    D_last = D;
  end  % of "while" for "run another case"
  disp('Model has been entered in SS form')
  model_entry = -1;                    % force model_entry menu to reappear
end  % of "if model_entry ==1"

%------------------ activate keyboard -----------------------
if model_entry == 2 ,
  disp('keyboard is active to allow:')
  disp('  * examining workspace')
  disp('  * entry of externally-defined model')
  disp('  * ... matrices must be in the arrays A B C D')
  disp('  * use of any Matlab commands')
  disp('To resume, enter "return"')
  keyboard
  model_entry = -1;                    % return to model_entry menu
end    % of "model_entry == 2"
%-------
end   % of " while model_entry....."
end   % of " if activity == 1"
%----------- end of model creation -------------------

%----------- Check controllability and observability -------------
if activity == 2,
  if isempty(A)~=1 | isempty(B)~=1     % check existance of matrices 
    disp(' ')
    disp('Check controllability:')
    [ns,ns] = size(A);
    Co = ctrb(A,B);
    c_eig = rank(Co);
    disp(' ')
    if c_eig == ns
	disp('The system is controllable.')
      else
        disp('The system is not controllable.')
	disp('The number of controllable eigenvalues is ')
        c_eig
    end
   else 
    disp('Cannot find state matrices A and/or B')
  end
  if isempty(A)~=1 | isempty(C)~=1   % check existance of matrices 
    disp(' ')
    disp('Check observability:')
    [ns,ns] = size(A);
    Ob = obsv(A,C);
    o_eig = rank(Ob);
    disp(' ')
    if o_eig == ns
	disp('The system is observable.')
      else
        disp('The system is not observable.')
	disp('The number of observable eigenvalues is ')
        o_eig
    end
   else 
    disp('Cannot find state matrices A and/or C')
  end
  activity = -1;                    % force activity menu to reappear
  disp('*****>');pause
end

%------------ end on check of controllability and observability ------------

%------------ Design full-state feedback gain using pole-placement ------------
if activity == 3,
  if isempty(A)~=1 | isempty(B)~=1  % check existance of matrices 
    disp(' ')
    disp('Design full-state feedback gain using pole-placement')
    p_last = [0.75 0.8 0.9]; 		% pole locations from example 10.2
    p_str = vec2str(p_last);
    pl = 'n';	% 'n' ==> pole location is not OK
    while pl == 'n',
      disp(' ')
      % allow user to enter new values for the desired poles
      p = input(['Enter poles "' p_str '"==> ']);
      if isempty(p), p = p_last; end   %  CR ==> use old value
      p_str = vec2str(p);  
      disp('Poles are: ')
      eval(['p = ' p_str]);
      pl = input('Poles OK? (y|n)','s');
      p_last = p;                    % save previous values in case want to reuse
    end   % of "while" for "run another case"
    disp(' ')
    F = place(A,B,p);
    disp('Full-state feedback gain is ')
    F 
   else 
    disp('Cannot find state matrices A and/or B')
  end
  disp('*****>');pause
  activity = -1;                     % force activity menu to reappear
  disp(' ')
end    % of "if activity == 3"
%------------------------------------------------------------

%------------ Design observer gain using pole-placement ---------
if activity == 4,
  if isempty(A)~=1 | isempty(C)~=1   % check existance of matrices 
    disp(' ')
    disp('Design observer gain using pole-placement')
    po_last = [0.5 0.6 0.75]; 		 % pole locations from example 10.5
    po_str = vec2str(po_last);
    pl = 'n';	% 'n' ==> pole location is not OK
    while pl == 'n',
      disp(' ')
      % allow user to enter new values for the desired poles
      po = input(['Enter poles "' po_str '"==> ']);
      if isempty(po), po = po_last; end   %  CR ==> use old value
      po_str = vec2str(po);  
      disp('Poles are: ')
      eval(['po = ' po_str]);
      pl = input('Poles OK? (y|n)','s');
      po_last = po;                       % save previous values in case want to reuse
    end   % of "while" for "run another case"
    disp(' ')
    L = place(A',C',po)';
    disp('Observer gain is ')
    L 
   else 
    disp('Cannot find state matrices A and/or C')
  end
  disp('*****>');pause
  activity = -1;   % force activity menu to reappear
  disp(' ')
end    % of "if activity == 4"
%------------------------------------------------

%------ Form observer controller and see its frequency response -------
%------             Also form closed-loop system                -------
if activity == 5,
  if isempty(A)~=1 | isempty(B)~=1 | isempty(C)~=1 | isempty(F)~=1 | ...
     isempty(L)~=1   % check existance of matrices 
    disp(' ')
    disp('Using F and L from previous design')
    disp(' ')
    disp('Controller system matrix')
    A_f = A-B*F-L*C                        % controller system matrix
    disp('*****>'), pause
    disp('controller poles')
    f_poles = eig(A_f)                     % controller poles
    disp('*****>'), pause
    disp('controller zeros')
    f_zeros = tzero(A_f,-L,-F,0)           % controller zeros
    disp('*****>'), pause    
    disp('Making controller Bode plot')
    Goc = ss(A_f,-L,-F,0,Ts)
    bode(Goc)                               % frequency response plot
    title('Frequency response of controller') 
    keyboard
    G = ss(A,B,C,D,Ts)
    G_cl = feedback(G*Goc,1)
    disp('Closed-loop system poles') 
    cl_loop_poles = pole(G_cl)             % closed-loop system poles
    disp('*****>');pause
    clf
    lfg = dcgain(G_cl);
    N = 1/lfg;
    G_clN = G_cl*N;
   else 
    disp('Cannot find state matrices A, B, C, F and/or L')
    subplot(1,1,1)
  end
  activity = -1;   % force activity menu to reappear
  disp(' ')
end    % of "if activity == 5"
%------------------------------------------------


%------------ Initial-condition response ------------
if activity == 6,
  if isempty(G_clN)~=1 % check existance of matrices 
    disp(' ')
    disp('Plot initial-condition response')
    z0_last = [1 -0.75 0.5 0 0 0]'
    z0_str = vec2str(z0_last);
    ic = 'n';	% 'n' ==> initial condition is not OK
    while ic == 'n',
      disp(' ')
      % allow user to enter new values for the initial condition
      z0 = input(['Enter initial condition "' z0_str '"==> ']);
      if isempty(z0), z0 = z0_last; end   %  CR ==> use old value
      z0_str = vec2str(z0);  
      disp('Initial condition is: ')
      eval(['z0 = ' z0_str]);
      ic = input('Initial condition OK? (y|n)','s');
      z0_last = z0;    % save previous values in case want to reuse
    end   % of "while" for "run another case"
    disp(' ')
    tf_last = 4;
    tf_str = vec2str(tf_last);
    tfin = 'n';	% 'n' ==> final time is not OK
    while tfin == 'n',
      disp(' ')
      % allow user to enter new values for the final time
      tf = input(['Enter final time "' tf_str '"==> ']);
      if isempty(tf), tf = tf_last; end   %  CR ==> use old value
      tf_str = vec2str(tf);  
      disp('final time is: ')
      eval(['tf = ' tf_str]);
      tfin = input('final time OK? (y|n)','s');
      tf_last = tf;                       % save previous values in case want to reuse
    end   % of "while" for "run another case"
    disp(' ')
    t = [0:Ts:tf]';
    u = 0*t;
    [y,t,x] = lsim(G_clN,u,t,z0);
    ns = length(z0);
    plot(t,x(:,1:ns/2),t,x(:,ns/2+1:ns),'--'), grid
    title('Initial condition response')
    disp('Solid lines are observer states, dashed lines ...')
    disp('  are system states')
    keyboard
   else 
    disp('Cannot find closed-loop systme G_clN')
  end
  activity = -1;                          % force activity menu to reappear
  disp(' ')
end    % of "if activity == 6"
%------------------------------------------------

%----------------- Step response ----------------
if activity == 7,
  if isempty(G_clN)~=1                    % check existance of matrices 
    title('Step input response')
    step(G_clN);
    disp(' ')
    disp('Compute time-domain perfermance measures')
    [ys,t,x] = step(G_clN);
    [Mo,tp,tr,ts,ess] = tstats(t,ys,1.0)
    keyboard
   else
    disp('Cannot find closed-loop system GclN')
  end
  activity = -1;                          % force activity menu to reappear
  disp(' ')
end    % of "if activity == 7"
%------------------------------------

end  % of "while" for activity menu
%%%%%%%%%%