CODEWITHSUNDEEP
matlabbisectionnorm
Dustin Malcolm
Dustin Malcolm

Reputation: 31

Problems with MATLAB nested statements and bisection

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

Upvotes: 3

Views: 115

Answers (1)

qbzenker
qbzenker

Reputation: 4652

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;

Upvotes: 0

Related Questions