function muBurlCE52(AutoFlag)

% mu synthesis for Burl computer exercise 5.2
% INPUTS: AutoFlag = true for automatic mu synthesis

if ~exist('AutoFlag', 'var')
    AutoFlag = false;
end

Integrator = tf(1, [1 0]) * eye(4);
sigmaxinv = tf(0.5, [1 0.1]) * eye(2);
A = [-10 0 0 0; 0 -10 0 0; 6 -3 -5 0; 0 0.5 0 -5];
B = [10 0; 0 10; 0 0; 0 0];
Bd = [0 0 0; 0 0 0; 1 1 0; 0 0 1];
C = [0 0 1 0; 0 0 0 1];
Cd = [2 0 0 0; 0 1 0 0; 0 0.4 0 0];

systemnames = 'Integrator A B Bd C Cd sigmaxinv';
inputvar = '[udel(3); w(2); u(2)]';
outputvar = '[Cd; sigmaxinv; w-C]';
input_to_Integrator = '[A+B+Bd]';
input_to_A = '[Integrator]';
input_to_B = '[u]';
input_to_Bd = '[udel]';
input_to_C = '[Integrator]';
input_to_Cd = '[Integrator]';
input_to_sigmaxinv = '[C-w]';
sysoutname = 'P';
sysic; % P = transfer function matrix
P = minreal(ss(P)); % P = state space model
% initialize mu synthesis variables
omega = logspace(-2, 2);
blk = [1 1; 1 1; 1 1; 2 2]; % 3x3 diagonal uncertainty, plus 2x2 performance block
nmeas = 2; % dimension of v
nu = 2; % dimension of u

if AutoFlag
    Delta = diag([ultidyn('D1',1), ultidyn('D2',1), ultidyn('D3',1)]);
    Punc = lft(Delta, P);
    opt = dkitopt('FrequencyVector', omega);
    [K, Clp, mubound, dkinfo] = dksyn(Punc, nmeas, nu, opt);
    disp([num2str(length(dkinfo)), ' iterations']);
    dkinfo = dkinfo(end);
    dkinfo = cell2mat(dkinfo);
    K = dkinfo.K;
    disp(['controller order = ', num2str(length(K.a))]);
    disp(['mu bound = ', num2str(dkinfo.Bnd)]);
    return;
end

DL = tf(eye(size(P))); % initial D scales
DR = DL;
while true
    % Find the H-inf controller with given D scales
    [K, N, gamma, info] = hinfsyn(DR*P*inv(DL), nmeas, nu, 'display', 'on');
    if isempty(K), break, end
    % Compute the SSV bounds of the system
    N = lft(P, K); % closed loop system (state space object)
    N = frd(N, omega); % frequency response data model of N
    [SSV, info] = mussv(N, blk); % SSV bounds of N
    mu = norm(SSV(1,1), inf); % upper bound of SSV of N
    ans = input(['mu = ', num2str(mu), ', controller order = ', num2str(length(K.a)), ' - another iteration (y/n)? '], 's');
    if ans ~= 'y', break, end;
    [DLfreq, DRfreq] = mussvunwrap(info); % fit D scales to SSV bounds
    for i = 1 : length(DL)-nu
        DLstable = genphase(DLfreq(i,i)); % stable min-phase freq response with magnitude DL(i,i)
        DRstable = genphase(DRfreq(i,i)); % stable min-phase freq response with magnitude DR(i,i)
        n = 3; % order of approximation to use
        DLss = fitfrd(DLstable, n); % nth order state space model that fits DL(i,i)
        DRss = fitfrd(DRstable, n); % nth order state space model that fits DR(i,i)
        DL(i,i) = tf(DLss); % convert DLss to transfer function form
        DR(i,i) = tf(DRss); % convert DRss to transfer function form
    end
end