% 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.