Mercurial > hg > octave-lojdl
changeset 3440:c487fb85b843
[project @ 2000-01-14 22:00:19 by jwe]
author | jwe |
---|---|
date | Fri, 14 Jan 2000 22:00:20 +0000 |
parents | 3234a698073a |
children | 36ae9880c594 |
files | scripts/control/base/dgkfdemo.m scripts/control/base/is_dgkf.m scripts/control/hinf/dgkfdemo.m scripts/control/hinf/is_dgkf.m |
diffstat | 4 files changed, 612 insertions(+), 612 deletions(-) [+] |
line wrap: on
line diff
deleted file mode 100644 --- a/scripts/control/base/dgkfdemo.m +++ /dev/null @@ -1,354 +0,0 @@ -## Copyright (C) 1996 Auburn University. All rights reserved. -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by the -## Free Software Foundation; either version 2, or (at your option) any -## later version. -## -## Octave is distributed in the hope that it will be useful, but WITHOUT -## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -## for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, write to the Free -## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {} dgkfdemo () -## Octave Controls toolbox demo: H2/Hinfinity options demos -## @end deftypefn - -## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> -## Created: June 1995 - -function dgkfdemo () - - save_val = page_screen_output; - page_screen_output = 1; - while (1) - clc - sel = 0; - while (sel > 10 || sel < 1) - sel = menu ("Octave H2/Hinfinity options demo", - "LQ regulator", - "LG state estimator", - "LQG optimal control design", - "H2 gain of a system", - "H2 optimal controller of a system", - "Hinf gain of a system", - "Hinf optimal controller of a SISO system", - "Hinf optimal controller of a MIMO system", - "Discrete-time Hinf optimal control by bilinear transform", - "Return to main demo menu"); - endwhile - if (sel == 1) - disp("Linear/Quadratic regulator design:") - disp("Compute optimal state feedback via the lqr command...") - help lqr - disp(" ") - disp("Example:") - A = [0, 1; -2, -1] - B = [0; 1] - Q = [1, 0; 0, 0] - R = 1 - disp("Q = state penalty matrix; R = input penalty matrix") - prompt - disp("Compute state feedback gain k, ARE solution P, and closed-loop") - disp("poles as follows:"); - cmd = "[k, p, e] = lqr(A,B,Q,R)"; - run_cmd - prompt - disp("A similar approach can be used for LTI discrete-time systems") - disp("by using the dlqr command in place of lqr (see LQG example).") - elseif (sel == 2) - disp("Linear/Gaussian estimator design:") - disp("Compute optimal state estimator via the lqe command...") - help lqe - disp(" ") - disp("Example:") - A = [0, 1; -2, -1] - disp("disturbance entry matrix G") - G = eye(2) - disp("Output measurement matrix C") - C = [0, 1] - SigW = [1, 0; 0, 1] - SigV = 1 - disp("SigW = input disturbance intensity matrix;") - disp("SigV = measurement noise intensity matrix") - prompt - disp("Compute estimator feedback gain k, ARE solution P, and estimator") - disp("poles via the command: ") - cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)"; - run_cmd - disp("A similar approach can be used for LTI discrete-time systems") - disp("by using the dlqe command in place of lqe (see LQG example).") - elseif (sel == 3) - disp("LQG optimal controller of a system:") - disp("Input accepted as either A,B,C matrices or in system data structure form") - disp("in both discrete and continuous time.") - disp("Example 1: continuous time design:") - prompt - help lqg - disp("Example system") - A = [0, 1; .5, .5]; - B = [0; 2]; - G = eye(2) - C = [1, 1]; - sys = ss2sys(A, [B, G], C); - sys = syssetsignals(sys,"in", ... - ["control input"; "disturbance 1"; "disturbance 2"]); - sysout(sys) - prompt - disp("Filtering/estimator parameters:") - SigW = eye(2) - SigV = 1 - prompt - disp("State space (LQR) parameters Q and R are:") - Q = eye(2) - R = 1 - cmd = "[K,Q1,P1,Ee,Er] = lqg(sys,SigW,SigV,Q,R,1);"; - run_cmd - disp("Check: closed loop system A-matrix is") - disp(" [A, B*Cc]") - disp(" [Bc*C, Ac ]") - cmd = "[Ac, Bc, Cc] = sys2ss(K);"; - run_cmd - cmd = "Acl = [A, B*Cc; Bc*C, Ac]"; - run_cmd - disp("Check: poles of Acl:") - Acl_poles = sortcom(eig(Acl)) - disp("Predicted poles from design = union(Er,Ee)") - cmd = "pred_poles = sortcom([Er; Ee])"; - run_cmd - disp("Example 2: discrete-time example") - cmd1 = "Dsys = ss2sys(A, [G, B], C, [0, 0, 0], 1);"; - cmd2 = "[K,Q1,P1,Ee,Er] = lqg(Dsys,SigW, SigV,Q,R);"; - disp("Run commands:") - cmd = cmd1; - run_cmd - cmd = cmd2; - run_cmd - prompt - disp("Check: closed loop system A-matrix is") - disp(" [A, B*Cc]") - disp(" [Bc*C, Ac ]") - [Ac,Bc,Cc] = sys2ss(K); - Acl = [A, B*Cc; Bc*C, Ac] - prompt - disp("Check: poles of Acl:") - Acl_poles = sortcom(eig(Acl)) - disp("Predicted poles from design = union(Er,Ee)") - pred_poles = sortcom([Er;Ee]) - elseif (sel == 4) - disp("H2 gain of a system: (Energy in impulse response)") - disp("Example 1: Stable plant:") - cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; - run_cmd - disp("Put into Packed system form:") - cmd = "Asys = ss2sys(A,B,C);"; - run_cmd - disp("Evaluate system 2-norm (impulse response energy):"); - cmd = "AsysH2 = h2norm(Asys)"; - run_cmd - disp("Compare with a plot of the system impulse response:") - tt = 0:0.1:20; - for ii=1:length(tt) - ht(ii) = C*expm(A*tt(ii))*B; - endfor - plot(tt,ht) - title("impulse response of example plant") - prompt - disp("Example 2: unstable plant") - cmd = "A = [0, 1; 2, 1]"; - eval(cmd); - cmd = "B = [0; 1]"; - eval(cmd); - cmd = "C = [1, 0]"; - eval(cmd); - cmd = "sys_poles = eig(A)"; - run_cmd - prompt - disp("Put into system data structure form:") - cmd="Bsys = ss2sys(A,B,C);"; - run_cmd - disp("Evaluate 2-norm:") - cmd = "BsysH2 = h2norm(Bsys)"; - run_cmd - disp(" ") - prompt("NOTICE: program returns a value without an error signal.") - disp("") - - elseif (sel == 5) - disp("H2 optimal controller of a system: command = h2syn:") - prompt - help h2syn - prompt - disp("Example system: double integrator with output noise and") - disp("input disturbance:") - disp(" "); - disp(" -------------------->y2"); - disp(" | _________"); - disp("u(t)-->o-->| 1/s^2 |-->o-> y1"); - disp(" ^ --------- ^"); - disp(" | |"); - disp(" w1(t) w2(t)"); - disp(" ") - disp("w enters the system through B1, u through B2") - disp("z = [y1; y2] is obtained through C1, y=y1 through C2"); - disp(" ") - cmd = "A = [0, 1; 0, 0]; B1 = [0, 0; 1, 0]; B2 = [0; 1];"; - disp(cmd) - eval(cmd); - cmd = "C1 = [1, 0; 0, 0]; C2 = [1, 0]; D11 = zeros(2);"; - disp(cmd) - eval(cmd); - cmd = "D12 = [0; 1]; D21 = [0, 1]; D22 = 0; D = [D11, D12; D21, D22];"; - disp(cmd) - eval(cmd); - disp("Design objective: compute U(s)=K(s)Y1(s) to minimize the closed") - disp("loop impulse response from w(t) =[w1; w2] to z(t) = [y1; y2]"); - prompt - disp("First: pack system:") - cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);"; - run_cmd - disp("Open loop multivariable Bode plot: (will take a moment)") - cmd="bode(Asys);"; - run_cmd - prompt("Press a key to close plot and continue"); - closeplot - disp("Controller design command: (only need 1st two output arguments)") - cmd="[K,gain, Kc, Kf, Pc, Pf] = h2syn(Asys,1,1);"; - run_cmd - disp("Controller is:") - cmd = "sysout(K)"; - run_cmd - disp(["returned gain value is: ",num2str(gain)]); - disp("Check: close the loop and then compute h2norm:") - prompt - cmd="K_loop = sysgroup(Asys,K);"; - run_cmd - cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);"; - run_cmd - cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);"; - run_cmd - cmd="gain_Kcl = h2norm(Kcl)"; - run_cmd - cmd="gain_err = gain_Kcl - gain"; - run_cmd - disp("Check: multivarible bode plot:") - cmd="bode(Kcl);"; - run_cmd - prompt - disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") - disp(" is_observable, is_detectable") - elseif (sel == 6) - disp("Hinfinity gain of a system: (max gain over all j-omega)") - disp("Example 1: Stable plant:") - cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; - run_cmd - disp("Pack into system format:") - cmd = "Asys = ss2sys(A,B,C);"; - run_cmd - disp("The infinity norm must be computed iteratively by") - disp("binary search. For this example, we select tolerance tol = 0.01, ") - disp("min gain gmin = 1e-2, max gain gmax=1e4.") - disp("Search quits when upper bound <= (1+tol)*lower bound.") - cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;"; - run_cmd - cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)" - run_cmd - disp("Check: look at max value of magntude Bode plot of Asys:"); - [M,P,w] = bode(Asys); - xlabel("Omega") - ylabel("|Asys(j omega)| ") - grid(); - semilogx(w,M); - disp(["Max magnitude is ",num2str(max(M)), ... - ", compared with gmin=",num2str(gmin)," and gmax=", ... - num2str(gmax),"."]) - prompt - disp("Example 2: unstable plant") - cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; - run_cmd - disp("Pack into system format:") - cmd = "Bsys = ss2sys(A,B,C);"; - run_cmd - disp("Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)") - BsysH2 = hinfnorm(Bsys,tol,gmin,gmax) - disp(" ") - disp("NOTICE: program returns a value without an error signal.") - disp("") - - elseif (sel == 7) - disp("Hinfinity optimal controller of a system: command = hinfsyn:") - prompt - help hinfsyn - prompt - disp("Example system: double integrator with output noise and") - disp("input disturbance:") - A = [0, 1; 0, 0] - B1 = [0, 0; 1, 0] - B2 = [0; 1] - C1 = [1, 0; 0, 0] - C2 = [1, 0] - D11 = zeros(2); - D12 = [0; 1]; - D21 = [0, 1]; - D22 = 0; - D = [D11, D12; D21, D22] - prompt - disp("First: pack system:") - cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);"; - run_cmd - prompt - disp("Open loop multivariable Bode plot: (will take a moment)") - cmd="bode(Asys);"; - run_cmd - prompt - disp("Controller design command: (only need 1st two output arguments)") - gmax = 1000 - gmin = 0.1 - gtol = 0.01 - cmd="[K,gain] = hinfsyn(Asys,1,1,gmin,gmax,gtol);"; - run_cmd - disp("Check: close the loop and then compute h2norm:") - prompt - cmd="K_loop = sysgroup(Asys,K);"; - run_cmd - cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);"; - run_cmd - cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);"; - run_cmd - cmd="gain_Kcl = hinfnorm(Kcl)"; - run_cmd - cmd="gain_err = gain_Kcl - gain"; - run_cmd - disp("Check: multivarible bode plot:") - cmd="bode(Kcl);"; - run_cmd - prompt - disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") - disp(" is_observable, is_detectable, buildssic") - elseif (sel == 8) - disp("Hinfinity optimal controller of MIMO system: command = hinfsyn:") - prompt - help hinfsyn - prompt - disp("Example system: Boeing 707-321 airspeed/pitch angle control") - disp(" ") - hinfdemo - elseif (sel == 9) - disp("Discrete time H-infinity control via bilinear transform"); - prompt - dhinfdemo - elseif (sel == 10) - return - endif - prompt - endwhile - page_screen_output = save_val; - -endfunction
deleted file mode 100644 --- a/scripts/control/base/is_dgkf.m +++ /dev/null @@ -1,258 +0,0 @@ -## Copyright (C) 1996, 1998 Auburn University. All rights reserved. -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by the -## Free Software Foundation; either version 2, or (at your option) any -## later version. -## -## Octave is distributed in the hope that it will be useful, but WITHOUT -## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -## for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, write to the Free -## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{retval}, @var{dgkf_struct} ] =} is_dgkf (@var{Asys}, @var{nu}, @var{ny}, @var{tol} ) -## Determine whether a continuous time state space system meets -## assumptions of DGKF algorithm. -## Partitions system into: -## @example -## [dx/dt] = [A | Bw Bu ][w] -## [ z ] [Cz | Dzw Dzu ][u] -## [ y ] [Cy | Dyw Dyu ] -## @end example -## or similar discrete-time system. -## If necessary, orthogonal transformations @var{Qw}, @var{Qz} and nonsingular -## transformations @var{Ru}, @var{Ry} are applied to respective vectors -## @var{w}, @var{z}, @var{u}, @var{y} in order to satisfy DGKF assumptions. -## Loop shifting is used if @var{Dyu} block is nonzero. -## -## @strong{Inputs} -## @table @var -## @item Asys -## system data structure -## @item nu -## number of controlled inputs -## @item ny -## number of measured outputs -## @item tol -## threshhold for 0. Default: 200@var{eps} -## @end table -## @strong{Outputs} -## @table @var -## @item retval -## true(1) if system passes check, false(0) otherwise -## @item dgkf_struct -## data structure of @code{is_dgkf} results. Entries: -## @table @var -## @item nw -## @itemx nz -## dimensions of @var{w}, @var{z} -## @item A -## system @var{A} matrix -## @item Bw -## (@var{n} x @var{nw}) @var{Qw}-transformed disturbance input matrix -## @item Bu -## (@var{n} x @var{nu}) @var{Ru}-transformed controlled input matrix; -## -## @strong{Note} @math{B = [Bw Bu] } -## @item Cz -## (@var{nz} x @var{n}) Qz-transformed error output matrix -## @item Cy -## (@var{ny} x @var{n}) @var{Ry}-transformed measured output matrix -## -## @strong{Note} @math{C = [Cz; Cy] } -## @item Dzu -## @item Dyw -## off-diagonal blocks of transformed @var{D} matrix that enter -## @var{z}, @var{y} from @var{u}, @var{w} respectively -## @item Ru -## controlled input transformation matrix -## @item Ry -## observed output transformation matrix -## @item Dyu_nz -## nonzero if the @var{Dyu} block is nonzero. -## @item Dyu -## untransformed @var{Dyu} block -## @item dflg -## nonzero if the system is discrete-time -## @end table -## @end table -## @code{is_dgkf} exits with an error if the system is mixed -## discrete/continuous -## -## @strong{References} -## @table @strong -## @item [1] -## Doyle, Glover, Khargonekar, Francis, "State Space Solutions -## to Standard H2 and Hinf Control Problems," IEEE TAC August 1989 -## @item [2] -## Maciejowksi, J.M.: "Multivariable feedback design," -## @end table -## @end deftypefn - -## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> -## Updated by John Ingram July 1996 to accept structured systems - -## Revised by Kai P. Mueller April 1998 to solve the general H_infinity -## problem using unitary transformations Q (on w and z) -## and non-singular transformations R (on u and y) such -## that the Dzu and Dyw matrices of the transformed plant -## -## ~ -## P (the variable Asys here) -## -## become -## -## ~ -1 T -## D = Q D R = [ 0 I ] or [ I ], -## 12 12 12 12 -## -## ~ T -## D = R D Q = [ 0 I ] or [ I ]. -## 21 21 21 21 -## -## This transformation together with the algorithm in [1] solves -## the general problem (see [2] for example). - -function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol) - - if (nargin < 3) | (nargin > 4) - usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})"); - elseif (! is_scalar(nu) | ! is_scalar(ny) ) - error("is_dgkf: arguments 2 and 3 must be scalars") - elseif (! is_struct(Asys) ) - error("Argument 1 must be a system data structure"); - endif - if(nargin < 4) - tol = 200*eps; - elseif( !is_sample(tol) ) - error("is_dgkf: tol must be a positive scalar") - endif - - retval = 1; # assume passes test - - dflg = is_digital(Asys); - [Anc, Anz, nin, nout ] = sysdimensions(Asys); - - if( Anz == 0 & Anc == 0 ) - error("is_dgkf: no system states"); - elseif( nu >= nin ) - error("is_dgkf: insufficient number of disturbance inputs"); - elseif( ny >= nout ) - error("is_dgkf: insufficient number of regulated outputs"); - endif - - nw = nin - nu; nw1 = nw + 1; - nz = nout - ny; nz1 = nz + 1; - - [A,B,C,D] = sys2ss(Asys); - ## scale input/output for numerical reasons - if(norm (C, "fro") * norm (B, "fro") == 0) - error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs"); - endif - xx = sqrt(norm(B, Inf) / norm(C, Inf)); - B = B / xx; C = C * xx; - - ## partition matrices - Bw = B(:,1:nw); Bu = B(:,nw1:nin); - Cz = C(1:nz,:); Dzw = D(1:nz,1:nw); Dzu = D(1:nz,nw1:nin); - Cy = C(nz1:nout,:); Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin); - - ## Check for loopo shifting - Dyu_nz = (norm(Dyu,Inf) != 0); - if (Dyu_nz) - warning("is_dgkf: D22 nonzero; performing loop shifting"); - endif - - ## 12 - rank condition at w = 0 - xx =[A, Bu; Cz, Dzu]; - [nr, nc] = size(xx); - irank = rank(xx); - if (irank != nc) - retval = 0; - warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ... - irank,nc,(Anc+Anz),nz,nu)); - warning(" *** 12-rank condition violated at w = 0."); - endif - - ## 21 - rank condition at w = 0 - xx =[A, Bw; Cy, Dyw]; - [nr, nc] = size(xx); - irank = rank(xx); - if (irank != nr) - retval = 0; - warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ... - irank,nr,(Anc+Anz),ny,nw)); - warning(" *** 21-rank condition violated at w = 0."); - endif - - ## can Dzu be transformed to become [0 I]' or [I]? - ## This ensures a normalized weight - [Qz, Ru] = qr(Dzu); - irank = rank(Ru); - if (irank != nu) - retval = 0; - warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank)); - warning(" *** Dzu does not have full column rank."); - endif - if (nu >= nz) - Qz = Qz(:,1:nu)'; - else - Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]'; - endif - Ru = Ru(1:nu,:); - - ## can Dyw be transformed to become [0 I] or [I]? - ## This ensures a normalized weight - [Qw, Ry] = qr(Dyw'); - irank = rank(Ry); - if (irank != ny) - retval = 0; - warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank)); - warning(" *** Dyw does not have full row rank."); - endif - - if (ny >= nw) - Qw = Qw(:,1:ny); - else - Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)]; - endif - Ry = Ry(1:ny,:)'; - - ## transform P by Qz/Ru and Qw/Ry - Bw = Bw*Qw; - Bu = Bu/Ru; - B = [Bw, Bu]; - Cz = Qz*Cz; - Cy = Ry\Cy; - C = [Cz; Cy]; - Dzw = Qz*Dzw*Qw; - Dzu = Qz*Dzu/Ru; - Dyw = Ry\Dyw*Qw; - - ## pack the return structure - dgkf_struct.nw = nw; - dgkf_struct.nu = nu; - dgkf_struct.nz = nz; - dgkf_struct.ny = ny; - dgkf_struct.A = A; - dgkf_struct.Bw = Bw; - dgkf_struct.Bu = Bu; - dgkf_struct.Cz = Cz; - dgkf_struct.Cy = Cy; - dgkf_struct.Dzw = Dzw; - dgkf_struct.Dzu = Dzu; - dgkf_struct.Dyw = Dyw; - dgkf_struct.Dyu = Dyu; - dgkf_struct.Ru = Ru; - dgkf_struct.Ry = Ry; - dgkf_struct.Dyu_nz = Dyu_nz; - dgkf_struct.dflg = dflg; - -endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/hinf/dgkfdemo.m @@ -0,0 +1,354 @@ +## Copyright (C) 1996 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by the +## Free Software Foundation; either version 2, or (at your option) any +## later version. +## +## Octave is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +## for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} dgkfdemo () +## Octave Controls toolbox demo: H2/Hinfinity options demos +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: June 1995 + +function dgkfdemo () + + save_val = page_screen_output; + page_screen_output = 1; + while (1) + clc + sel = 0; + while (sel > 10 || sel < 1) + sel = menu ("Octave H2/Hinfinity options demo", + "LQ regulator", + "LG state estimator", + "LQG optimal control design", + "H2 gain of a system", + "H2 optimal controller of a system", + "Hinf gain of a system", + "Hinf optimal controller of a SISO system", + "Hinf optimal controller of a MIMO system", + "Discrete-time Hinf optimal control by bilinear transform", + "Return to main demo menu"); + endwhile + if (sel == 1) + disp("Linear/Quadratic regulator design:") + disp("Compute optimal state feedback via the lqr command...") + help lqr + disp(" ") + disp("Example:") + A = [0, 1; -2, -1] + B = [0; 1] + Q = [1, 0; 0, 0] + R = 1 + disp("Q = state penalty matrix; R = input penalty matrix") + prompt + disp("Compute state feedback gain k, ARE solution P, and closed-loop") + disp("poles as follows:"); + cmd = "[k, p, e] = lqr(A,B,Q,R)"; + run_cmd + prompt + disp("A similar approach can be used for LTI discrete-time systems") + disp("by using the dlqr command in place of lqr (see LQG example).") + elseif (sel == 2) + disp("Linear/Gaussian estimator design:") + disp("Compute optimal state estimator via the lqe command...") + help lqe + disp(" ") + disp("Example:") + A = [0, 1; -2, -1] + disp("disturbance entry matrix G") + G = eye(2) + disp("Output measurement matrix C") + C = [0, 1] + SigW = [1, 0; 0, 1] + SigV = 1 + disp("SigW = input disturbance intensity matrix;") + disp("SigV = measurement noise intensity matrix") + prompt + disp("Compute estimator feedback gain k, ARE solution P, and estimator") + disp("poles via the command: ") + cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)"; + run_cmd + disp("A similar approach can be used for LTI discrete-time systems") + disp("by using the dlqe command in place of lqe (see LQG example).") + elseif (sel == 3) + disp("LQG optimal controller of a system:") + disp("Input accepted as either A,B,C matrices or in system data structure form") + disp("in both discrete and continuous time.") + disp("Example 1: continuous time design:") + prompt + help lqg + disp("Example system") + A = [0, 1; .5, .5]; + B = [0; 2]; + G = eye(2) + C = [1, 1]; + sys = ss2sys(A, [B, G], C); + sys = syssetsignals(sys,"in", ... + ["control input"; "disturbance 1"; "disturbance 2"]); + sysout(sys) + prompt + disp("Filtering/estimator parameters:") + SigW = eye(2) + SigV = 1 + prompt + disp("State space (LQR) parameters Q and R are:") + Q = eye(2) + R = 1 + cmd = "[K,Q1,P1,Ee,Er] = lqg(sys,SigW,SigV,Q,R,1);"; + run_cmd + disp("Check: closed loop system A-matrix is") + disp(" [A, B*Cc]") + disp(" [Bc*C, Ac ]") + cmd = "[Ac, Bc, Cc] = sys2ss(K);"; + run_cmd + cmd = "Acl = [A, B*Cc; Bc*C, Ac]"; + run_cmd + disp("Check: poles of Acl:") + Acl_poles = sortcom(eig(Acl)) + disp("Predicted poles from design = union(Er,Ee)") + cmd = "pred_poles = sortcom([Er; Ee])"; + run_cmd + disp("Example 2: discrete-time example") + cmd1 = "Dsys = ss2sys(A, [G, B], C, [0, 0, 0], 1);"; + cmd2 = "[K,Q1,P1,Ee,Er] = lqg(Dsys,SigW, SigV,Q,R);"; + disp("Run commands:") + cmd = cmd1; + run_cmd + cmd = cmd2; + run_cmd + prompt + disp("Check: closed loop system A-matrix is") + disp(" [A, B*Cc]") + disp(" [Bc*C, Ac ]") + [Ac,Bc,Cc] = sys2ss(K); + Acl = [A, B*Cc; Bc*C, Ac] + prompt + disp("Check: poles of Acl:") + Acl_poles = sortcom(eig(Acl)) + disp("Predicted poles from design = union(Er,Ee)") + pred_poles = sortcom([Er;Ee]) + elseif (sel == 4) + disp("H2 gain of a system: (Energy in impulse response)") + disp("Example 1: Stable plant:") + cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; + run_cmd + disp("Put into Packed system form:") + cmd = "Asys = ss2sys(A,B,C);"; + run_cmd + disp("Evaluate system 2-norm (impulse response energy):"); + cmd = "AsysH2 = h2norm(Asys)"; + run_cmd + disp("Compare with a plot of the system impulse response:") + tt = 0:0.1:20; + for ii=1:length(tt) + ht(ii) = C*expm(A*tt(ii))*B; + endfor + plot(tt,ht) + title("impulse response of example plant") + prompt + disp("Example 2: unstable plant") + cmd = "A = [0, 1; 2, 1]"; + eval(cmd); + cmd = "B = [0; 1]"; + eval(cmd); + cmd = "C = [1, 0]"; + eval(cmd); + cmd = "sys_poles = eig(A)"; + run_cmd + prompt + disp("Put into system data structure form:") + cmd="Bsys = ss2sys(A,B,C);"; + run_cmd + disp("Evaluate 2-norm:") + cmd = "BsysH2 = h2norm(Bsys)"; + run_cmd + disp(" ") + prompt("NOTICE: program returns a value without an error signal.") + disp("") + + elseif (sel == 5) + disp("H2 optimal controller of a system: command = h2syn:") + prompt + help h2syn + prompt + disp("Example system: double integrator with output noise and") + disp("input disturbance:") + disp(" "); + disp(" -------------------->y2"); + disp(" | _________"); + disp("u(t)-->o-->| 1/s^2 |-->o-> y1"); + disp(" ^ --------- ^"); + disp(" | |"); + disp(" w1(t) w2(t)"); + disp(" ") + disp("w enters the system through B1, u through B2") + disp("z = [y1; y2] is obtained through C1, y=y1 through C2"); + disp(" ") + cmd = "A = [0, 1; 0, 0]; B1 = [0, 0; 1, 0]; B2 = [0; 1];"; + disp(cmd) + eval(cmd); + cmd = "C1 = [1, 0; 0, 0]; C2 = [1, 0]; D11 = zeros(2);"; + disp(cmd) + eval(cmd); + cmd = "D12 = [0; 1]; D21 = [0, 1]; D22 = 0; D = [D11, D12; D21, D22];"; + disp(cmd) + eval(cmd); + disp("Design objective: compute U(s)=K(s)Y1(s) to minimize the closed") + disp("loop impulse response from w(t) =[w1; w2] to z(t) = [y1; y2]"); + prompt + disp("First: pack system:") + cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);"; + run_cmd + disp("Open loop multivariable Bode plot: (will take a moment)") + cmd="bode(Asys);"; + run_cmd + prompt("Press a key to close plot and continue"); + closeplot + disp("Controller design command: (only need 1st two output arguments)") + cmd="[K,gain, Kc, Kf, Pc, Pf] = h2syn(Asys,1,1);"; + run_cmd + disp("Controller is:") + cmd = "sysout(K)"; + run_cmd + disp(["returned gain value is: ",num2str(gain)]); + disp("Check: close the loop and then compute h2norm:") + prompt + cmd="K_loop = sysgroup(Asys,K);"; + run_cmd + cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);"; + run_cmd + cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);"; + run_cmd + cmd="gain_Kcl = h2norm(Kcl)"; + run_cmd + cmd="gain_err = gain_Kcl - gain"; + run_cmd + disp("Check: multivarible bode plot:") + cmd="bode(Kcl);"; + run_cmd + prompt + disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") + disp(" is_observable, is_detectable") + elseif (sel == 6) + disp("Hinfinity gain of a system: (max gain over all j-omega)") + disp("Example 1: Stable plant:") + cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; + run_cmd + disp("Pack into system format:") + cmd = "Asys = ss2sys(A,B,C);"; + run_cmd + disp("The infinity norm must be computed iteratively by") + disp("binary search. For this example, we select tolerance tol = 0.01, ") + disp("min gain gmin = 1e-2, max gain gmax=1e4.") + disp("Search quits when upper bound <= (1+tol)*lower bound.") + cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;"; + run_cmd + cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)" + run_cmd + disp("Check: look at max value of magntude Bode plot of Asys:"); + [M,P,w] = bode(Asys); + xlabel("Omega") + ylabel("|Asys(j omega)| ") + grid(); + semilogx(w,M); + disp(["Max magnitude is ",num2str(max(M)), ... + ", compared with gmin=",num2str(gmin)," and gmax=", ... + num2str(gmax),"."]) + prompt + disp("Example 2: unstable plant") + cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)"; + run_cmd + disp("Pack into system format:") + cmd = "Bsys = ss2sys(A,B,C);"; + run_cmd + disp("Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)") + BsysH2 = hinfnorm(Bsys,tol,gmin,gmax) + disp(" ") + disp("NOTICE: program returns a value without an error signal.") + disp("") + + elseif (sel == 7) + disp("Hinfinity optimal controller of a system: command = hinfsyn:") + prompt + help hinfsyn + prompt + disp("Example system: double integrator with output noise and") + disp("input disturbance:") + A = [0, 1; 0, 0] + B1 = [0, 0; 1, 0] + B2 = [0; 1] + C1 = [1, 0; 0, 0] + C2 = [1, 0] + D11 = zeros(2); + D12 = [0; 1]; + D21 = [0, 1]; + D22 = 0; + D = [D11, D12; D21, D22] + prompt + disp("First: pack system:") + cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);"; + run_cmd + prompt + disp("Open loop multivariable Bode plot: (will take a moment)") + cmd="bode(Asys);"; + run_cmd + prompt + disp("Controller design command: (only need 1st two output arguments)") + gmax = 1000 + gmin = 0.1 + gtol = 0.01 + cmd="[K,gain] = hinfsyn(Asys,1,1,gmin,gmax,gtol);"; + run_cmd + disp("Check: close the loop and then compute h2norm:") + prompt + cmd="K_loop = sysgroup(Asys,K);"; + run_cmd + cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);"; + run_cmd + cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);"; + run_cmd + cmd="gain_Kcl = hinfnorm(Kcl)"; + run_cmd + cmd="gain_err = gain_Kcl - gain"; + run_cmd + disp("Check: multivarible bode plot:") + cmd="bode(Kcl);"; + run_cmd + prompt + disp("Related functions: is_dgkf, is_controllable, is_stabilizable,") + disp(" is_observable, is_detectable, buildssic") + elseif (sel == 8) + disp("Hinfinity optimal controller of MIMO system: command = hinfsyn:") + prompt + help hinfsyn + prompt + disp("Example system: Boeing 707-321 airspeed/pitch angle control") + disp(" ") + hinfdemo + elseif (sel == 9) + disp("Discrete time H-infinity control via bilinear transform"); + prompt + dhinfdemo + elseif (sel == 10) + return + endif + prompt + endwhile + page_screen_output = save_val; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/hinf/is_dgkf.m @@ -0,0 +1,258 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by the +## Free Software Foundation; either version 2, or (at your option) any +## later version. +## +## Octave is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +## for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{retval}, @var{dgkf_struct} ] =} is_dgkf (@var{Asys}, @var{nu}, @var{ny}, @var{tol} ) +## Determine whether a continuous time state space system meets +## assumptions of DGKF algorithm. +## Partitions system into: +## @example +## [dx/dt] = [A | Bw Bu ][w] +## [ z ] [Cz | Dzw Dzu ][u] +## [ y ] [Cy | Dyw Dyu ] +## @end example +## or similar discrete-time system. +## If necessary, orthogonal transformations @var{Qw}, @var{Qz} and nonsingular +## transformations @var{Ru}, @var{Ry} are applied to respective vectors +## @var{w}, @var{z}, @var{u}, @var{y} in order to satisfy DGKF assumptions. +## Loop shifting is used if @var{Dyu} block is nonzero. +## +## @strong{Inputs} +## @table @var +## @item Asys +## system data structure +## @item nu +## number of controlled inputs +## @item ny +## number of measured outputs +## @item tol +## threshhold for 0. Default: 200@var{eps} +## @end table +## @strong{Outputs} +## @table @var +## @item retval +## true(1) if system passes check, false(0) otherwise +## @item dgkf_struct +## data structure of @code{is_dgkf} results. Entries: +## @table @var +## @item nw +## @itemx nz +## dimensions of @var{w}, @var{z} +## @item A +## system @var{A} matrix +## @item Bw +## (@var{n} x @var{nw}) @var{Qw}-transformed disturbance input matrix +## @item Bu +## (@var{n} x @var{nu}) @var{Ru}-transformed controlled input matrix; +## +## @strong{Note} @math{B = [Bw Bu] } +## @item Cz +## (@var{nz} x @var{n}) Qz-transformed error output matrix +## @item Cy +## (@var{ny} x @var{n}) @var{Ry}-transformed measured output matrix +## +## @strong{Note} @math{C = [Cz; Cy] } +## @item Dzu +## @item Dyw +## off-diagonal blocks of transformed @var{D} matrix that enter +## @var{z}, @var{y} from @var{u}, @var{w} respectively +## @item Ru +## controlled input transformation matrix +## @item Ry +## observed output transformation matrix +## @item Dyu_nz +## nonzero if the @var{Dyu} block is nonzero. +## @item Dyu +## untransformed @var{Dyu} block +## @item dflg +## nonzero if the system is discrete-time +## @end table +## @end table +## @code{is_dgkf} exits with an error if the system is mixed +## discrete/continuous +## +## @strong{References} +## @table @strong +## @item [1] +## Doyle, Glover, Khargonekar, Francis, "State Space Solutions +## to Standard H2 and Hinf Control Problems," IEEE TAC August 1989 +## @item [2] +## Maciejowksi, J.M.: "Multivariable feedback design," +## @end table +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Updated by John Ingram July 1996 to accept structured systems + +## Revised by Kai P. Mueller April 1998 to solve the general H_infinity +## problem using unitary transformations Q (on w and z) +## and non-singular transformations R (on u and y) such +## that the Dzu and Dyw matrices of the transformed plant +## +## ~ +## P (the variable Asys here) +## +## become +## +## ~ -1 T +## D = Q D R = [ 0 I ] or [ I ], +## 12 12 12 12 +## +## ~ T +## D = R D Q = [ 0 I ] or [ I ]. +## 21 21 21 21 +## +## This transformation together with the algorithm in [1] solves +## the general problem (see [2] for example). + +function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol) + + if (nargin < 3) | (nargin > 4) + usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})"); + elseif (! is_scalar(nu) | ! is_scalar(ny) ) + error("is_dgkf: arguments 2 and 3 must be scalars") + elseif (! is_struct(Asys) ) + error("Argument 1 must be a system data structure"); + endif + if(nargin < 4) + tol = 200*eps; + elseif( !is_sample(tol) ) + error("is_dgkf: tol must be a positive scalar") + endif + + retval = 1; # assume passes test + + dflg = is_digital(Asys); + [Anc, Anz, nin, nout ] = sysdimensions(Asys); + + if( Anz == 0 & Anc == 0 ) + error("is_dgkf: no system states"); + elseif( nu >= nin ) + error("is_dgkf: insufficient number of disturbance inputs"); + elseif( ny >= nout ) + error("is_dgkf: insufficient number of regulated outputs"); + endif + + nw = nin - nu; nw1 = nw + 1; + nz = nout - ny; nz1 = nz + 1; + + [A,B,C,D] = sys2ss(Asys); + ## scale input/output for numerical reasons + if(norm (C, "fro") * norm (B, "fro") == 0) + error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs"); + endif + xx = sqrt(norm(B, Inf) / norm(C, Inf)); + B = B / xx; C = C * xx; + + ## partition matrices + Bw = B(:,1:nw); Bu = B(:,nw1:nin); + Cz = C(1:nz,:); Dzw = D(1:nz,1:nw); Dzu = D(1:nz,nw1:nin); + Cy = C(nz1:nout,:); Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin); + + ## Check for loopo shifting + Dyu_nz = (norm(Dyu,Inf) != 0); + if (Dyu_nz) + warning("is_dgkf: D22 nonzero; performing loop shifting"); + endif + + ## 12 - rank condition at w = 0 + xx =[A, Bu; Cz, Dzu]; + [nr, nc] = size(xx); + irank = rank(xx); + if (irank != nc) + retval = 0; + warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ... + irank,nc,(Anc+Anz),nz,nu)); + warning(" *** 12-rank condition violated at w = 0."); + endif + + ## 21 - rank condition at w = 0 + xx =[A, Bw; Cy, Dyw]; + [nr, nc] = size(xx); + irank = rank(xx); + if (irank != nr) + retval = 0; + warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ... + irank,nr,(Anc+Anz),ny,nw)); + warning(" *** 21-rank condition violated at w = 0."); + endif + + ## can Dzu be transformed to become [0 I]' or [I]? + ## This ensures a normalized weight + [Qz, Ru] = qr(Dzu); + irank = rank(Ru); + if (irank != nu) + retval = 0; + warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank)); + warning(" *** Dzu does not have full column rank."); + endif + if (nu >= nz) + Qz = Qz(:,1:nu)'; + else + Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]'; + endif + Ru = Ru(1:nu,:); + + ## can Dyw be transformed to become [0 I] or [I]? + ## This ensures a normalized weight + [Qw, Ry] = qr(Dyw'); + irank = rank(Ry); + if (irank != ny) + retval = 0; + warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank)); + warning(" *** Dyw does not have full row rank."); + endif + + if (ny >= nw) + Qw = Qw(:,1:ny); + else + Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)]; + endif + Ry = Ry(1:ny,:)'; + + ## transform P by Qz/Ru and Qw/Ry + Bw = Bw*Qw; + Bu = Bu/Ru; + B = [Bw, Bu]; + Cz = Qz*Cz; + Cy = Ry\Cy; + C = [Cz; Cy]; + Dzw = Qz*Dzw*Qw; + Dzu = Qz*Dzu/Ru; + Dyw = Ry\Dyw*Qw; + + ## pack the return structure + dgkf_struct.nw = nw; + dgkf_struct.nu = nu; + dgkf_struct.nz = nz; + dgkf_struct.ny = ny; + dgkf_struct.A = A; + dgkf_struct.Bw = Bw; + dgkf_struct.Bu = Bu; + dgkf_struct.Cz = Cz; + dgkf_struct.Cy = Cy; + dgkf_struct.Dzw = Dzw; + dgkf_struct.Dzu = Dzu; + dgkf_struct.Dyw = Dyw; + dgkf_struct.Dyu = Dyu; + dgkf_struct.Ru = Ru; + dgkf_struct.Ry = Ry; + dgkf_struct.Dyu_nz = Dyu_nz; + dgkf_struct.dflg = dflg; + +endfunction