3

New to the site, but I'm trying to hone up on some MATLAB skills for work and school, and was looking for some help with the following:

I want to write my own algorithm for finding the Hinf norm of a system using bisection, like the MATLAB function 'hinfsyn' does.

I've included the code I have so far:

function [ hnorm ] = matmath(A,B,C,D,glow,ghigh,tol)

if 2*(ghigh-glow) < tol
    gam  = (ghigh+glow)/2;
    hnorm = gam;
else
    Dgam = ((gam^2)*eye(size(D)))-(D'*D);
    A_clp = [ A + (B/Dgam*D'*C) -B/Dgam*B'
        (C'*C) + C'*D/Dgam*D'*C -A'-(C'*D/Dgam*B')];
    eigcl = eig(A_clp);
    for i = 1:length(eigcl)
        if real(eig(i)) == 0
            glow = gam;
        else
            ghigh = gam;
        end
    end
end

I've rationalized the problem to a few steps:

  1. With gamma bounds used as an input, compute the first iteration: gam = (ghigh-glow)/2. If 2*(ghigh-glow) < tol, then the program stops with hinf = gam.
  2. Compute eignevalues of closed loop A matrix.
  3. Check for purely imaginary eigenvalues. If there exists a purely imaginary eigenvalue, new glow = gam. Otherwise, set ghigh = gam.
  4. Continue iterating until there gamma tolerance is satisfied.

I believe my matrix calculations are correct, but I'm having a hard time with the if/for statements. My code only completes the first iteration when I execute it in MATLAB. Any tips would be greatly appreciated! Thank you!

Update: here is the code I've simplified, which successfully completes one iteration and has been tested for many different values.

function [ hnorm ] = matmath(A,B,C,D,glow,ghigh,tol)

gam  = (ghigh+glow)/2;

if 2*(ghigh-glow) < tol
    hnorm = gam
else    

Dgam = ((gam^2)*eye(size(D)))-(D'*D);
A_clp = [ A + (B/Dgam*D'*C) -B/Dgam*B'
        (C'*C) + C'*D/Dgam*D'*C -A'-(C'*D/Dgam*B')];
eig_clp = eig(A_clp)

for z = 1:length(eig_clp)
    if abs(real(eig_clp(z)))<1e-10
        glow = gam
        break
    end
end

ghigh = gam

end
Dustin Malcolm
  • 31
  • 3

1 Answers1

0

If you want your code to run until you achieve 2*(gmax - gmin) < tol, use a while loop:

function [ hnorm ] = matmath(A,B,C,D,glow,ghigh,tol)

while 2*(ghigh-glow) >= tol
    Dgam = ((gam^2)*eye(size(D)))-(D'*D);
    A_clp = [ A + (B/Dgam*D'*C) -B/Dgam*B'
        (C'*C) + C'*D/Dgam*D'*C -A'-(C'*D/Dgam*B')];
    eigcl = eig(A_clp);
    for i = 1:length(eigcl)
        if real(eig(i)) == 0
            glow = gam;
        else
            ghigh = gam;
        end
    end
end

gam  = (ghigh+glow)/2;
hnorm = gam;
qbzenker
  • 4,412
  • 3
  • 16
  • 23
  • Thank you for your help! That form makes more sense. However, my code still stops after the first iteration. I have a suspicion that the for loop checking for purely imaginary eigenvalues is preventing gam from being updated. – Dustin Malcolm May 01 '17 at 19:35
  • So, I realized that my input arguments for glow and ghigh were flipped, which is why my code was stopping on the first iteration. However, when I flipped the responses, i would get an error for "Undefined function or variable "gam". I moved that statement to the top of the program, and then it ran. Gamma then becomes updated at the for loop, but stops after the reassignment of gam, and doesn't return to the top – Dustin Malcolm May 01 '17 at 20:02
  • OK, my guess is that something about the calculations is messed up. – qbzenker May 01 '17 at 20:08
  • I see an error in my code. "if real(eig(i)) == 0" SHOULD BE "if real(eigcl(i)) == 0." However, I also believe that's where my code is getting hung up. I've run example eigenvalues through that statement, and even ones that have a real value of 1e-12 won't meet the 0.0 requirement, even though that precision meets what I desire in that statement. Is there a better way to check that inequality? I assume round-off with floating point numbers is causing an issue. – Dustin Malcolm May 01 '17 at 21:28
  • definitely could be - you could try something like `if abs(real(eigcl(i))) < 1e-8` – qbzenker May 01 '17 at 21:31
  • So the abs() and < 1e-10 seemed to work. However, when I run the program now, MATLAB gets stuck in a loop and won't stop running off of the values from the first iteration. I've simplified my code and removed that while statement to verify that ghigh or glow are activated on the appropriate eigenvalues for the first iteration, but I just need to figure out a good way for the new ghigh or glow to get passed back to the top for the next iteration. I'll edit my original question above to include what I have so far. – Dustin Malcolm May 01 '17 at 22:21