% ice_root.m

% Adapt ice_root.m to do calculations for your problem sets.

% ice_root calculates the RGA (relative gain array) and DC (disturbance
% cost) for the 10.492 heat exchanger network

% Run ice_root from the Matlab command line.  Its use requires that a 
% column vector called "xguess" be available in memory, and that 3 files
% be available in the user's directory:

% reduced_Newton.m - presented by Prof. Beers in 10.34, it applies the
% Newton-Raphson method (with some embellishments) to find the root of a
% system of nonlinear algebraic equations.  It has been modified for 
% 10.492, and should need no further attention from the user.

% calc_f.m - calculates the system of equations, given a guess of the root.
% When calc_f.m returns a vector of zeroes, the root has been found. 
% calc_f.m is used by reduced_Newton.m.  The user should modify the 
% calc_f.m template to contain the network equations.  Any parameters
% that are to be held constant (such as Cp and U) may be set there.

% calc_Jac.m - calculates the Jacobian matrix (multidimensional slope) of
% the equation system in calc_f.m.  It should need no modification by the
% user.

disp ('**********  Computations for the 10.492 HX Network  *************')

% set up the rootfinding parameters in a data structure
% their significance is described in the reduced_Newton.m file
Options.max_iter = 200;
Options.max_iter_LS = 200;
Options.rtol = 1;
Options.atol = 1.e-5;
Options.step_tol = 1.e-3;
Options.verbose = 1;
Options.use_range = 0;
Options.range = 0;

disp ('The control parameters for root-finding routine reduced_Newton.m')
disp (Options)

% set up the numerical differentiation spacing
% dy/dx ~ (y((1+xdiff)*x) - y((1-xdiff)*x))/(2*x*xdiff)
% vary this value in repeated runs to check whether the calculated 
% derivatives are sensitive to the value chosen for xdiff
xdiff = 0.01;

disp ('The differentiation interval as a percentage of the ref value')
disp (xdiff)

%determine whether the user should form the column vector that contains the
%initial guess

if exist('xguess') == 0
    disp ('You must first define the initial guess column vector "xguess"')
    disp ('The length of xguess must be consistent with calc_f.m')
    return
end

% determine whether the structure Param exists.
% if so, use it, because the user may wish to modify it from the command
% line.
% if not, create it here
if exist('Param') == 0
    Param.w1 = 8.2;
    Param.w2 = 9.1;
    Param.w3 = 5.6;
    Param.T1 = 250;
    Param.T2 = 136;
    Param.T3 = 100;
end

% set the manipulated variable reference conditions
w2r = Param.w2;
w3r = Param.w3;
disp ('Manipulated variable reference values')
w2r, w3r

% set the disturbance variable reference conditions
w1r = Param.w1;
T1r = Param.T1;
T2r = Param.T2;
T3r = Param.T3;
disp ('Disturbance variable reference values')
w1r, T1r, T2r, T3r

% prepare linear system matrices
% set the matrix dimensions according to the numbers of variables
Pm = zeros(3,2);  % 3 CV; 2 MV
Pd = zeros(3,4);  % 3 CV; 4 DV

% prepare scaling matrices for the variables
% insert appropriate scaling ranges as required for individual variables
S = 2*blkdiag(10, 10, 10);      % all temperatures at +/- 10 deg C
Sm = 2*blkdiag(w2r,  w3r);      % all flows at +/- the reference value
Sd = 2*blkdiag(w1r,  10,  10,  10);

% prepare a summary matrix to hold DC calculations
% its length is number-of-values ^ number-of-dist-variables
% its width is 1 + number-of-dist-variables + number-of-manip-variables
DC_summary = zeros(81,7);

% Section 1: calculate the reference condition CV from the 
% reference DV and MV
% ----------------------------------------------------------------
% call reduced_Newton to find the root satisfying the variables of Param
[xref,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xguess,@calc_f,@calc_Jac,Options,Param); 

% display output
disp ('Calculation of the Reference Case')
disp ('number of iterations'); disp(iter_conv)

% if this reference case is satisfactorily converged, the other
% calculations, which are conducted nearby, are likely to be successful.
% for normal convergence, clean up imaginary components that may have
% arisen during root-finding.
% leave them to attract attention otherwise.
if iflag == 1
    disp ('convergence normal')
    xref = real(xref);
else 
    disp ('did not converge; run with improved guess')
end

% xref is the reference case root; ensure that the definition is
% consistent between this program and calc_f.m
disp ('controlled variables T6r T8r T9r'); disp(xref(1:3))
disp ('intermediate variables T4r T5r T7r'); disp (xref(4:6))


% Section 2: calculate matrix Pm by successive perturbation of the
% manipulated variables.  Each differentiation calculates one column of Pm.
% -----------------------------------------------------------
% w/r/t MV W2
Param.w2 = w2r*(1+xdiff);  %positive perturbation of w2
[rootplus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Param.w2 = w2r*(1-xdiff);  %negative perturbation of w2
[rootminus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Pm(1:3,1) = (rootplus(1:3)-rootminus(1:3))/(2*xdiff*w2r);
Param.w2 = w2r;  %return w2 to reference value

% w/r/t MV W3
Param.w3 = w3r*(1+xdiff);
[rootplus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Param.w3 = w3r*(1-xdiff);
[rootminus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Pm(1:3,2) = (rootplus(1:3)-rootminus(1:3))/(2*xdiff*w3r);
Param.w3 = w3r;

disp (' ')
disp ('The matrix Pm')
disp (Pm)

% Section 3: calculate matrix Pd by successive perturbation of the
% disturbance variables.  Each differentiation calculates one column of Pd.
% -----------------------------------------------------------
% w/r/t DV W1
Param.w1 = w1r*(1+xdiff);
[rootplus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Param.w1 = w1r*(1-xdiff);
[rootminus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Pd(1:3,1) = (rootplus(1:3)-rootminus(1:3))/(2*xdiff*w1r);
Param.w1 = w1r;

% w/r/t DV T1
Param.T1 = T1r*(1+xdiff);
[rootplus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Param.T1 = T1r*(1-xdiff);
[rootminus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Pd(1:3,2) = (rootplus(1:3)-rootminus(1:3))/(2*xdiff*T1r);
Param.T1 = T1r;

% w/r/t DV T2
Param.T2 = T2r*(1+xdiff);
[rootplus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Param.T2 = T2r*(1-xdiff);
[rootminus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Pd(1:3,3) = (rootplus(1:3)-rootminus(1:3))/(2*xdiff*T2r);
Param.T2 = T2r;

% w/r/t DV T3
Param.T3 = T3r*(1+xdiff);
[rootplus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Param.T3 = T3r*(1-xdiff);
[rootminus,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 

Pd(1:3,4) = (rootplus(1:3)-rootminus(1:3))/(2*xdiff*T3r);
Param.T3 = T3r;

disp (' ')
disp ('The matrix Pd')
disp (Pd)

% Section 4: scale the coefficient matrices
% ---------------------------------------------------------
disp (' ')
disp ('The scaled coefficient matrices')
Pmstar = S\Pm*Sm
Pdstar = S\Pd*Sd

% Section 5: calculate RGA 
% because we have only 2 MV for 3 CV, we must calculate several
% combinations.
% separate cases would be unnecessary for balanced CV and MV
% -----------------------------------------------------------
% case 1: control T6 and T8
disp ('control T6 and T8')
Pmcase = Pm(1:2,:)      % select appropriate rows
RGA = Pmcase.*inv(Pmcase)'

% case 2: control T6 and T9
disp ('control T6 and T9')
Pmcase = [Pm(1,:);Pm(3,:)]      % select appropriate rows
RGA = Pmcase.*inv(Pmcase)'

% case 3: control T8 and T9
disp ('control T8 and T9')
Pmcase = Pm(2:3,:)      % select appropriate rows
RGA = Pmcase.*inv(Pmcase)'

% Section 6: calculate DC
% case 3 has been selected for calculating DC; this, of course, required
% coding this section after preliminary runs had been made of the RGA.
% ----------------------------------------------------------------
% select the subset of Pmstar and Pdstar to use for DC
Pmsel = Pmstar(2:3,:);
Pdsel = Pdstar(2:3,:);

% calculate the quotient of Pdsel and Pmsel
Pquo = Pmsel\Pdsel;

% set disturbance limits (add variables as required)
xd1 = [ -0.25  0  0.25 ]; % variable w1
xd2 = [  -0.5  0  0.5 ]; % variable T1
xd3 = [  -0.5  0  0.5 ]; % variable T2
xd4 = [  -0.5  0  0.5 ]; % variable T3

% loop to calculate combinations of the DV (add loops as required)
case_num = 0;
for i1 = 1:3                %perturb variable xd1
    for i2 = 1:3            %perturb variable xd2
        for i3 = 1:3        %perturb variable xd3
            for i4 = 1:3    %perturb variable xd4
                case_num = case_num+1;
                xdstar = [xd1(i1); xd2(i2); xd3(i3); xd4(i4)];
                xmstar = -Pquo*xdstar;
                DC = norm(xmstar);
                DC_summary (case_num,:) = [DC xdstar' xmstar'];
            end
        end
    end
end
% find the maximum DC
DC_sort = sortrows(DC_summary);

disp ('largest values of Disturbance Cost')
disp ('     DC   Dist_Variables                             Manip_Variables')
length = size(DC_sort,1);
disp (DC_sort ((length-5):length,:))

% Section 7: at the maximum DC
% examine the required manipulated variable capacity, and compare
% the linear and nonlinear models at the extreme disturbance values
% ----------------------------------------------------------------
% compose a vector of disturbance variables at the worst DC
xdstar_max = DC_sort(length,2:5)';

% compose a vector of manipulated variables at this condition
xmstar_pc = DC_sort(length,6:7)';

% compute the controlled variables via the linear model
ystar_pc = Pmstar*xmstar_pc + Pdstar*xdstar_max;

% convert to deviation variables for output
y_pc = S*ystar_pc;
disp ('controlled variable deviation at worst disturbance - linear model (degC)')
disp (y_pc)

%convert the DV (max) and MV (pc) to deviation variables
xd_max = Sd*xdstar_max;
xm_pc = Sm*xmstar_pc;

% set the Param structure to the extreme conditions
Param.w1 = xd_max(1) + w1r;
Param.w2 = xm_pc(1) + w2r;
Param.w3 = xm_pc(2) + w3r;
Param.T1 = xd_max(2) + T1r;
Param.T2 = xd_max(3) + T2r;
Param.T3 = xd_max(4) + T3r;

% compute the controlled variables via the nonlinear model
[y_nonlin,iflag,iter_conv,x_traj,f_traj] = ... 
     reduced_Newton(xref,@calc_f,@calc_Jac,Options,Param); 
 
% express nonlinear output as deviation variables
y_pcnl = zeros(3,1);
y_pcnl(1) = y_nonlin(1) - xref(1);
y_pcnl(2) = y_nonlin(2) - xref(2);
y_pcnl(3) = y_nonlin(3) - xref(3);
disp ('controlled variable deviation at worst disturbance - nonlinear model (degC)')
disp (y_pcnl(1:3))

% return Param to reference values in case user wishes further runs
Param.w1 = w1r;
Param.w2 = w2r;
Param.w3 = w3r;
Param.T1 = T1r;
Param.T2 = T2r;
Param.T3 = T3r;

% express manipulated variable in physical form
xm_req(1) = xm_pc(1) + w2r;
xm_req(2) = xm_pc(2) + w3r;
disp ('check that manipulated variable motion is not extreme')
disp ('w2 (kg s-1)'); disp (xm_req(1))
disp ('w3 (kg s-1)'); disp (xm_req(2))

% Dr. Barry S. Johnston, Copyright 2004.