Mercurial > hg > octave-terminal
changeset 7131:a184bc985c37
[project @ 2007-11-08 15:55:02 by jwe]
author | jwe |
---|---|
date | Thu, 08 Nov 2007 15:55:02 +0000 |
parents | 5eeb46c784d7 |
children | b01db194c526 |
files | scripts/ChangeLog scripts/control/hinf/h2syn.m scripts/control/hinf/hinf_ctr.m scripts/control/hinf/hinfnorm.m scripts/control/hinf/hinfsyn.m scripts/control/hinf/hinfsyn_chk.m scripts/control/hinf/is_dgkf.m scripts/control/hinf/wgt1o.m |
diffstat | 8 files changed, 299 insertions(+), 239 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,10 @@ +2007-11-08 John W. Eaton <jwe@octave.org> + + * control/hinf/h2syn.m, control/hinf/hinf_ctr.m, + control/hinf/hinfnorm.m, control/hinf/hinfsyn.m, + control/hinf/hinfsyn_chk.m, control/hinf/is_dgkf.m, + control/hinf/wgt1o.m: Style fixes. + 2007-11-08 David Bateman <dbateman@free.fr> * plot/quiver.m: Fix arrowheads.
--- a/scripts/control/hinf/h2syn.m +++ b/scripts/control/hinf/h2syn.m @@ -79,43 +79,59 @@ function [K, gain, Kc, Kf, Pc, Pf] = h2syn (Asys, nu, ny, tol) - if ((nargin < 3) | (nargin > 4)) + if (nargin < 3 || nargin > 4) print_usage (); - elseif(nargin == 3 ) - [chkdgkf,dgs] = is_dgkf(Asys,nu,ny); - elseif(nargin == 4) - [chkdgkf,dgs] = is_dgkf(Asys,nu,ny,tol); + elseif (nargin == 3) + [chkdgkf, dgs] = is_dgkf (Asys, nu, ny); + elseif (nargin == 4) + [chkdgkf, dgs] = is_dgkf (Asys, nu, ny, tol); endif - if (!chkdgkf ) - disp("h2syn: system does not meet required assumptions") - help is_dgkf - error("h2syn: exit"); + if (! chkdgkf) + error ("h2syn: system does not meet required assumptions") endif ## extract dgs information - nw = dgs.nw; nu = dgs.nu; - A = dgs.A; Bw = dgs.Bw; Bu = dgs.Bu; - Cz = dgs.Cz; Dzw = dgs.Dzw; Dzu = dgs.Dzu; nz = dgs.nz; - Cy = dgs.Cy; Dyw = dgs.Dyw; Dyu = dgs.Dyu; ny = dgs.ny; + nw = dgs.nw; + nu = dgs.nu; + nz = dgs.nz; + ny = dgs.ny; + + A = dgs.A; + + Bw = dgs.Bw; + Bu = dgs.Bu; + + Cz = dgs.Cz; + Cy = dgs.Cy; + + Dzw = dgs.Dzw; + Dzu = dgs.Dzu; + + Dyw = dgs.Dyw; + Dyu = dgs.Dyu; + d22nz = dgs.Dyu_nz; + dflg = dgs.dflg; - if(norm(Dzw,Inf) > norm([Dzw, Dzu ; Dyw, Dyu],Inf)*1e-12) - warning("h2syn: Dzw nonzero; feedforward not implemented") + if (norm (Dzw, Inf) > norm ([Dzw, Dzu; Dyw, Dyu], Inf)*1e-12) + warning ("h2syn: Dzw nonzero; feedforward not implemented") Dzw D = [Dzw, Dzu ; Dyw, Dyu] endif ## recover i/o transformations - Ru = dgs.Ru; Ry = dgs.Ry; - [ncstates, ndstates, nout, nin] = sysdimensions(Asys); - Atsam = sysgettsam(Asys); - [Ast, Ain, Aout] = sysgetsignals(Asys); + Ru = dgs.Ru; + Ry = dgs.Ry; - if(dgs.dflg == 0) - Pc = are(A,Bu*Bu',Cz'*Cz); # solve control, filtering ARE's - Pf = are(A',Cy'*Cy,Bw*Bw'); + [ncstates, ndstates, nout, nin] = sysdimensions (Asys); + Atsam = sysgettsam (Asys); + [Ast, Ain, Aout] = sysgetsignals (Asys); + + if (dgs.dflg == 0) + Pc = are (A, Bu*Bu', Cz'*Cz); # solve control, filtering ARE's + Pf = are(A', Cy'*Cy, Bw*Bw'); F2 = -Bu'*Pc; # calculate feedback gains L2 = -Pf*Cy'; @@ -126,61 +142,62 @@ else ## discrete time solution - error("h2syn: discrete-time case not yet implemented") - Pc = dare(A,Bu*Bu',Cz'*Cz); - Pf = dare(A',Cy'*Cy,Bw*Bw'); + error ("h2syn: discrete-time case not yet implemented") + Pc = dare (A, Bu*Bu', Cz'*Cz); + Pf = dare (A', Cy'*Cy, Bw*Bw'); endif nn = ncstates + ndstates; - In = eye(nn); + In = eye (nn); KA = A + Bu*F2 + L2*Cy; - Kc1 = ss(AF2,Bw,CzF2,zeros(nz,nw)); - Kf1 = ss(AL2,BwL2,F2,zeros(nu,nw)); + Kc1 = ss (AF2, Bw, CzF2, zeros (nz, nw)); + Kf1 = ss (AL2, BwL2, F2, zeros (nu, nw)); - g1 = h2norm(Kc1); - g2 = h2norm(Kf1); + g1 = h2norm (Kc1); + g2 = h2norm (Kf1); ## compute optimal closed loop gain - gain = sqrt ( g1*g1 + g2*g2 ); + gain = sqrt (g1*g1 + g2*g2); - if(nargout) - Kst = strappend(Ast,"_K"); - Kin = strappend(Aout((nout-ny+1):(nout)),"_K"); - Kout = strappend(Ain((nin-nu+1):(nin)),"_K"); + if (nargout) + Kst = strappend (Ast, "_K"); + Kin = strappend (Aout((nout-ny+1):(nout)), "_K"); + Kout = strappend (Ain((nin-nu+1):(nin)), "_K"); ## compute systems for return - K = ss(KA,-L2/Ru,Ry\F2,zeros(nu,ny),Atsam,ncstates,ndstates,Kst,Kin,Kout); + K = ss (KA, -L2/Ru, Ry\F2, zeros(nu,ny), Atsam, ncstates, + ndstates, Kst, Kin, Kout); endif if (nargout > 2) ## system full information control state names - stname2 = strappend(Ast,"_FI"); + stname2 = strappend (Ast, "_FI"); ## system full information control input names - inname2 = strappend(Ast,"_FI_in"); + inname2 = strappend (Ast, "_FI_in"); ## system full information control output names - outname2 = strappend(Aout(1:(nout-ny)),"_FI_out"); + outname2 = strappend (Aout(1:(nout-ny)), "_FI_out"); nz = rows (Cz); nw = columns (Bw); - Kc = ss(AF2, In, CzF2, zeros(nz,nn), Atsam, ... - ncstates, ndstates, stname2, inname2, outname2); + Kc = ss (AF2, In, CzF2, zeros(nz,nn), Atsam, + ncstates, ndstates, stname2, inname2, outname2); endif if (nargout >3) ## fix system state estimator state names - stname3 = strappend(Ast,"_Kf"); + stname3 = strappend (Ast, "_Kf"); ## fix system state estimator input names - inname3 = strappend(Ast,"_Kf_noise"); + inname3 = strappend (Ast, "_Kf_noise"); ## fix system state estimator output names - outname3 = strappend(Ast,"_est"); + outname3 = strappend (Ast, "_est"); - Kf = ss(AL2, BwL2, In, zeros(nn,nw),Atsam, ... - ncstates, ndstates, stname3, inname3,outname3); + Kf = ss (AL2, BwL2, In, zeros(nn,nw),Atsam, + ncstates, ndstates, stname3, inname3, outname3); endif endfunction
--- a/scripts/control/hinf/hinf_ctr.m +++ b/scripts/control/hinf/hinf_ctr.m @@ -143,7 +143,7 @@ ## non-zero D22 is a special case if (d22nz) if (rank (eye(nu) + d11hat*D22) < nu) - error(" *** cannot compute controller for D22 non-zero."); + error (" *** cannot compute controller for D22 non-zero."); endif d22new = [D22, zeros(ny,ny); zeros(nu,nu), 0*D22'];
--- a/scripts/control/hinf/hinfnorm.m +++ b/scripts/control/hinf/hinfnorm.m @@ -115,93 +115,93 @@ function [g, gmin, gmax] = hinfnorm (sys, tol, gmin, gmax, ptol) - if((nargin == 0) || (nargin > 4)) + if (nargin == 0 || nargin > 4) print_usage (); - elseif(!isstruct(sys)) - error("Sys must be a system data structure"); + elseif (! isstruct (sys)) + error ("Sys must be a system data structure"); endif ## set defaults where applicable - if(nargin < 5) + if (nargin < 5) ptol = 1e-9; # pole tolerance endif - if(nargin < 4) + if (nargin < 4) gmax = 1e9; # max gain value endif - dflg = is_digital(sys); - sys = sysupdate(sys,"ss"); - [A,B,C,D] = sys2ss(sys); - [n,nz,m,p] = sysdimensions(sys); + dflg = is_digital (sys); + sys = sysupdate (sys, "ss"); + [A, B, C, D] = sys2ss (sys); + [n, nz, m, p] = sysdimensions (sys); ## eigenvalues of A must all be stable - if(!is_stable(sys)) - warning(["hinfnorm: unstable system (is_stable, ptol=",num2str(ptol), ... - "), returning Inf"]); + if (! is_stable (sys)) + warning ("hinfnorm: unstable system (is_stable, ptol=%g), returning Inf", + ptol); g = Inf; endif - Dnrm = norm(D); - if(nargin < 3) - gmin = max(1e-9,Dnrm); # min gain value - elseif(gmin < Dnrm) - warning(["hinfnorm: setting Gmin=||D||=",num2str(Dnrm)]); + Dnrm = norm (D); + if (nargin < 3) + gmin = max (1e-9, Dnrm); # min gain value + elseif (gmin < Dnrm) + warning ("hinfnorm: setting Gmin=||D||=%g", Dnrm); endif - if(nargin < 2) + if (nargin < 2) tol = 0.001; # convergence measure for gmin, gmax endif ## check for scalar input arguments 2...5 - if( ! (isscalar(tol) && isscalar(gmin) - && isscalar(gmax) && isscalar(ptol)) ) - error("hinfnorm: tol, gmin, gmax, ptol must be scalars"); + if (! isscalar (tol) && isscalar (gmin) + && isscalar (gmax) && isscalar (ptol)) ) + error ("hinfnorm: tol, gmin, gmax, ptol must be scalars"); endif - In = eye(n+nz); - Im = eye(m); - Ip = eye(p); + In = eye (n+nz); + Im = eye (m); + Ip = eye (p); ## find the Hinf norm via binary search - while((gmax/gmin - 1) > tol) + while (gmax/gmin - 1 > tol) g = (gmax+gmin)/2; - if(dflg) + if (dflg) ## multiply g's through in formulas to avoid extreme magnitudes... Rg = g^2*Im - D'*D; Ak = A + (B/Rg)*D'*C; Ck = g^2*C'*((g^2*Ip-D*D')\C); ## set up symplectic generalized eigenvalue problem per Iglesias & Glover - s1 = [Ak , zeros(nz) ; -Ck, In ]; - s2 = [In, -(B/Rg)*B' ; zeros(nz) , Ak' ]; + s1 = [Ak , zeros(nz); -Ck, In]; + s2 = [In, -(B/Rg)*B'; zeros(nz), Ak']; ## guard against roundoff again: zero out extremely small values ## prior to balancing s1 = s1 .* (abs(s1) > ptol*norm(s1,"inf")); s2 = s2 .* (abs(s2) > ptol*norm(s2,"inf")); - [cc,dd,s1,s2] = balance(s1,s2); - [qza,qzb,zz,pls] = qz(s1,s2,"S"); # ordered qz decomposition - eigerr = abs(abs(pls)-1); - normH = norm([s1,s2]); + [cc, dd, s1, s2] = balance (s1, s2); + [qza, qzb, zz, pls] = qz (s1, s2, "S"); # ordered qz decomposition + eigerr = abs (abs(pls)-1); + normH = norm ([s1, s2]); Hb = [s1, s2]; ## check R - B' X B condition (Iglesias and Glover's paper) X = zz((nz+1):(2*nz),1:nz)/zz(1:nz,1:nz); - dcondfailed = min(real( eig(Rg - B'*X*B)) < ptol); + dcondfailed = min (real (eig (Rg - B'*X*B)) < ptol); else Rinv = inv(g*g*Im - (D' * D)); - H = [A + B*Rinv*D'*C, B*Rinv*B'; ... + H = [A + B*Rinv*D'*C, B*Rinv*B'; -C'*(Ip + D*Rinv*D')*C, -(A + B*Rinv*D'*C)']; ## guard against roundoff: zero out extremely small values prior ## to balancing - H = H .* (abs(H) > ptol*norm(H,"inf")); - [DD,Hb] = balance(H); - pls = eig(Hb); - eigerr = abs(real(pls)); - normH = norm(H); + H = H .* (abs (H) > ptol * norm (H, "inf")); + [DD, Hb] = balance (H); + pls = eig (Hb); + eigerr = abs (real (pls)); + normH = norm (H); dcondfailed = 0; # digital condition; doesn't apply here endif - if( (min(eigerr) <= ptol * normH) | dcondfailed) + if ((min (eigerr) <= ptol * normH) | dcondfailed) gmin = g; else gmax = g;
--- a/scripts/control/hinf/hinfsyn.m +++ b/scripts/control/hinf/hinfsyn.m @@ -131,47 +131,64 @@ function [K, g, GW, Xinf, Yinf] = hinfsyn (Asys, nu, ny, gmin, gmax, gtol, ptol, tol) - if( (nargin < 1) | (nargin > 8) ) + if (nargin < 1 || nargin > 8) print_usage (); endif ## set default arguments - if(nargin < 8) + if (nargin < 8) tol = 200*eps; - elseif(!is_sample(tol)) - error("tol must be a positive scalar.") + elseif (! is_sample (tol)) + error ("tol must be a positive scalar.") endif - if(nargin < 7) + if (nargin < 7) ptol = 1e-9; - elseif(!is_sample(ptol)) - error("hinfsyn: ptol must be a positive scalar"); + elseif (! is_sample (ptol)) + error ("hinfsyn: ptol must be a positive scalar"); endif - if(!is_sample(gmax) | !is_sample(gmin) | !is_sample(gtol) ) - error(["hinfsyn: gmax=",num2str(gmax),", gmin=",num2str(gmin), ... - "gtol=",num2str(gtol), " must be positive scalars."]) + if (! is_sample (gmax) || ! is_sample (gmin) || ! is_sample (gtol)) + error ("hinfsyn: gmax=%g, gmin=%g, gtol=%g must be positive scalars", + gmax, gmin, gtol); endif - [chkdgkf,dgs] = is_dgkf(Asys,nu,ny,tol); + [chkdgkf, dgs] = is_dgkf (Asys, nu, ny, tol); if (! chkdgkf ) - disp("hinfsyn: system does not meet required assumptions") - help is_dgkf - error("hinfsyn: exit"); + error ("hinfsyn: system does not meet required assumptions") endif ## extract dgs information - nw = dgs.nw; nu = dgs.nu; - A = dgs.A; B1 = dgs.Bw; B2 = dgs.Bu; - C1 = dgs.Cz; D11 = dgs.Dzw; D12 = dgs.Dzu; nz = dgs.nz; - C2 = dgs.Cy; D21 = dgs.Dyw; D22 = dgs.Dyu; ny = dgs.ny; + nw = dgs.nw; + nu = dgs.nu; + nz = dgs.nz; + ny = dgs.ny; + + A = dgs.A; + + B1 = dgs.Bw; + B2 = dgs.Bu; + + C1 = dgs.Cz; + C2 = dgs.Cy; + + D11 = dgs.Dzw; + D12 = dgs.Dzu; + D21 = dgs.Dyw; + D22 = dgs.Dyu; + d22nz = dgs.Dyu_nz; + dflg = dgs.dflg; ## recover i/o transformations - R12 = dgs.Ru; R21 = dgs.Ry; - [ncstates, ndstates, nin, nout] = sysdimensions(Asys); - Atsam = sysgettsam(Asys); - [Ast, Ain, Aout] = sysgetsignals(Asys); + R12 = dgs.Ru; + R21 = dgs.Ry; + + [ncstates, ndstates, nin, nout] = sysdimensions (Asys); + + Atsam = sysgettsam (Asys); + + [Ast, Ain, Aout] = sysgetsignals (Asys); BB = [B1, B2]; CC = [C1 ; C2]; @@ -182,21 +199,19 @@ ## perform binary search to find gamma min ghi = gmax; ## derive a lower lower bound for gamma from D matrices - xx1 = norm((eye(nz) - (D12/(D12'*D12))*D12')*D11); - xx2 = norm(D11*(eye(nw)-(D21'/(D21*D21'))*D21)); - glo = max(xx1, xx2); + xx1 = norm ((eye(nz) - (D12/(D12'*D12))*D12')*D11); + xx2 = norm (D11*(eye(nw)-(D21'/(D21*D21'))*D21)); + glo = max (xx1, xx2); if (glo > gmin) - disp(" *** D matrices indicate a greater value of gamma min."); - fprintf(" gamma min (%f) superseeded by %f.", gmin, glo); + warning (" *** D matrices indicate a greater value of gamma min."); + warning (" gamma min (%f) superseeded by %f.", gmin, glo); glo = xx1; else glo = gmin; endif if (glo > ghi) - fprintf(" *** lower bound of gamma greater than upper bound(%f)", ... - glo, ghi); - disp(" *** unable to continue, Goodbye."); - return; + error (" *** lower bound of gamma greater than upper bound (%g, %g)", + glo, ghi); endif de = ghi - glo; @@ -242,118 +257,131 @@ ## now do the search - while (!iteration_finished) + while (! iteration_finished) switch (search_state) - case (0) g = ghi; - case (1) g = glo; - case (2) g = 0.5 * (ghi + glo); - otherwise error(" *** This should never happen!"); + case 0 + g = ghi; + case 1 + g = glo; + case 2 + g = 0.5 * (ghi + glo); + otherwise + error (" *** This should never happen!"); endswitch - printf("%10.4f ", g); + printf ("%10.4f ", g); ## computing R and R~ d1dot = [D11, D12]; - R = zeros(nin, nin); + R = zeros (nin, nin); R(1:nw,1:nw) = -g*g*eye(nw); R = R + d1dot' * d1dot; ddot1 = [D11; D21]; - Rtilde = zeros(nout, nout); + Rtilde = zeros (nout, nout); Rtilde(1:nz,1:nz) = -g*g*eye(nz); Rtilde = Rtilde + ddot1 * ddot1'; - [Xinf,x_ha_err] = hinfsyn_ric(A,BB,C1,d1dot,R,ptol); - [Yinf,y_ha_err] = hinfsyn_ric(A',CC',B1',ddot1',Rtilde,ptol); + [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol); + [Yinf, y_ha_err] = hinfsyn_ric (A', CC', B1', ddot1', Rtilde, ptol); ## assume failure for this gamma passed = 0; - rerr=""; - if (!x_ha_err && !y_ha_err) + rerr = ""; + if (! x_ha_err && ! y_ha_err) ## test spectral radius condition - rho = max(abs(eig(Xinf * Yinf))); + rho = max (abs (eig (Xinf * Yinf))); if (rho < g*g) ## spectral radius condition passed passed = 1; else - rerr = sprintf("rho=%f",rho); + rerr = sprintf ("rho=%f",rho); endif endif - if(x_ha_err >= 0 & x_ha_err <= 6) - printf("%s",errmesg{x_ha_err+1}); + if (x_ha_err >= 0 && x_ha_err <= 6) + printf ("%s", errmesg{x_ha_err+1}); xerr = errdesx{x_ha_err+1}; else - error(" *** Xinf fail: this should never happen!"); + error (" *** Xinf fail: this should never happen!"); endif - if(y_ha_err >= 0 & y_ha_err <= 6) - printf("%s",errmesg{y_ha_err+1}); + if (y_ha_err >= 0 && y_ha_err <= 6) + printf ("%s", errmesg{y_ha_err+1}); yerr = errdesy{y_ha_err+1}; else - error(" *** Yinf fail: this should never happen!"); + error (" *** Yinf fail: this should never happen!"); endif - if(passed) printf(" y all tests passed.\n"); - else printf(" n %s/%s%s\n",xerr,yerr,rerr); endif + if (passed) + printf (" y all tests passed.\n"); + else + printf (" n %s/%s%s\n", xerr, yerr, rerr); + endif if (passed && (de/g < gtol)) search_state = 3; endif switch (search_state) - case (0) - if (!passed) + case 0 + if (! passed) ## upper bound must pass but did not fprintf(" *** the upper bound of gamma (%f) is too small.\n", g); iteration_finished = 2; else search_state = 1; endif - case (1) - if (!passed) search_state = 2; + case 1 + if (! passed) + search_state = 2; else ## lower bound must not pass but passed - fprintf(" *** the lower bound of gamma (%f) passed.\n", g); + fprintf (" *** the lower bound of gamma (%f) passed.\n", g); iteration_finished = 3; endif - case (2) + case 2 ## Normal case; must check that singular R, Rtilde wasn't the problem. - if ((!passed) & (x_ha_err != 6) & (y_ha_err != 6) ) glo = g; - else ghi = g; endif + if (! passed && x_ha_err != 6 && y_ha_err != 6) + glo = g; + else + ghi = g; + endif de = ghi - glo; - case (3) iteration_finished = 1; # done - otherwise error(" *** This should never happen!"); + case 3 + iteration_finished = 1; # done + otherwise + error (" *** This should never happen!"); endswitch endwhile - printf("----------------------------------------"); - printf("--------------------------------------\n"); + printf ("----------------------------------------"); + printf ("--------------------------------------\n"); if (iteration_finished != 1) K = []; else ## success: compute controller - fprintf(" hinfsyn final: glo=%f ghi=%f, test gain g=%f\n", \ - glo, ghi, g); - printf("----------------------------------------"); - printf("--------------------------------------\n"); - Z = inv(eye(ncstates) - Yinf*Xinf/g/g); + fprintf (" hinfsyn final: glo=%f ghi=%f, test gain g=%f\n", + glo, ghi, g); + printf ("----------------------------------------"); + printf ("--------------------------------------\n"); + Z = inv (eye(ncstates) - Yinf*Xinf/g/g); F = -R \ (d1dot'*C1 + BB'*Xinf); H = -(B1*ddot1' + Yinf*CC') / Rtilde; - K = hinf_ctr(dgs,F,H,Z,g); + K = hinf_ctr (dgs, F, H, Z, g); - Kst = strappend(Ast,"_K"); - Kin = strappend(Aout((nout-ny+1):(nout)),"_K"); - Kout = strappend(Ain((nin-nu+1):(nin)),"_K"); - [Ac, Bc, Cc, Dc] = sys2ss(K); - K = ss(Ac,Bc,Cc,Dc,Atsam,ncstates,ndstates,Kst,Kin,Kout); + Kst = strappend (Ast, "_K"); + Kin = strappend (Aout((nout-ny+1):(nout)),"_K"); + Kout = strappend (Ain((nin-nu+1):(nin)),"_K"); + [Ac, Bc, Cc, Dc] = sys2ss (K); + K = ss (Ac, Bc, Cc, Dc, Atsam, ncstates, ndstates, Kst, Kin, Kout); if (nargout >= 3) - GW = starp(Asys, K); + GW = starp (Asys, K); endif endif - elseif(ndstates) + elseif (ndstates) ## discrete time solution - error("hinfsyn: discrete-time case not yet implemented") + error ("hinfsyn: discrete-time case not yet implemented") endif
--- a/scripts/control/hinf/hinfsyn_chk.m +++ b/scripts/control/hinf/hinfsyn_chk.m @@ -88,8 +88,8 @@ ## Construct the two Hamiltonians g2 = 1/(g*g); - Hc = [ A , g2*B1*B1' - B2*B2'; -C1'*C1 , -A']; - Hf = [ A' , g2*C1'*C1 - C2'*C2; -B1*B1' , -A]; + Hc = [A, g2*B1*B1' - B2*B2'; -C1'*C1, -A']; + Hf = [A', g2*C1'*C1 - C2'*C2; -B1*B1', -A]; ## check if Hc, Hf are in dom(Ric) Hcminval = min (abs (real (eig (Hc))));
--- a/scripts/control/hinf/is_dgkf.m +++ b/scripts/control/hinf/is_dgkf.m @@ -131,84 +131,95 @@ function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol) - if (nargin < 3) | (nargin > 4) + if (nargin < 3) || (nargin > 4) print_usage (); - elseif (! isscalar(nu) | ! isscalar(ny) ) - error("is_dgkf: arguments 2 and 3 must be scalars") - elseif (! isstruct(Asys) ) - error("Argument 1 must be a system data structure"); + elseif (! isscalar (nu) || ! isscalar (ny)) + error ("is_dgkf: arguments 2 and 3 must be scalars") + elseif (! isstruct (Asys)) + error ("Argument 1 must be a system data structure"); endif - if(nargin < 4) + if (nargin < 4) tol = 200*eps; - elseif( !is_sample(tol) ) - error("is_dgkf: tol must be a positive scalar") + 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); + 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"); + 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; + nw = nin - nu; + nz = nout - ny; - [A,B,C,D] = sys2ss(Asys); + nw1 = nw + 1; + 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"); + 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)); + 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); + Bw = B(:,1:nw); + Bu = B(:,nw1:nin); + + Cz = C(1:nz,:); + Cy = C(nz1:nout,:); + + Dzw = D(1:nz,1:nw); + Dzu = D(1:nz,nw1:nin); + + Dyw = D(nz1:nout,1:nw); + Dyu = D(nz1:nout,nw1:nin); ## Check for loopo shifting - Dyu_nz = (norm(Dyu,Inf) != 0); + Dyu_nz = (norm (Dyu, Inf) != 0); if (Dyu_nz) - warning("is_dgkf: D22 nonzero; performing loop shifting"); + 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); + [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."); + warning ("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); + [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."); + warning ("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); + [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."); + warning ("*** 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)'; @@ -219,12 +230,12 @@ ## can Dyw be transformed to become [0 I] or [I]? ## This ensures a normalized weight - [Qw, Ry] = qr(Dyw'); - irank = rank(Ry); + [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."); + warning ("*** rank(Dyw(%d x %d) = %d", ny, nw, irank); + warning (" *** Dyw does not have full row rank"); endif if (ny >= nw) @@ -246,22 +257,22 @@ 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; + 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
--- a/scripts/control/hinf/wgt1o.m +++ b/scripts/control/hinf/wgt1o.m @@ -59,20 +59,17 @@ print_usage (); endif - if(nargout > 1) - print_usage (); - endif - if (vl == vh) a = []; b = []; c = []; else - a = [-2*pi*fc]; - b = [-2*pi*fc]; - c = [vh-vl]; + a = -2*pi*fc; + b = -2*pi*fc; + c = vh-vl; endif - d=[vh]; + d = vh; - wsys = ss(a,b,c,d); + wsys = ss (a, b, c, d); + endfunction