% =========================================================================
%                        BROBEKK INCINERATION PLANT
% =========================================================================
% Nonlinear model of the Brobekk garbage incineration plant in Oslo solved 
% with the TOMLAB solver Minos.
% 
% The objective function is to minimize the cost related to pump and fan 
% work. Mathematically, this is formulated as
% 
%   J = min c'x             where c is a cost vector; and x the vector of
%                           optimizable variables (DOFs)
%
%   s.t.
%
%   A*x <= b                Linear nonequalities
%   Aeq*x = beq             Linear equalities
%   LB <= x <= UB           Lower and upper bounds
%   C(x) <= 0               Nonlinear inequalities
%   Ceq(x) = 0              Nonlinear equalities
%
% The degrees of freedom are:
%
%   x(1:5) = u          Valves
%   x(6:7) = P          Pump duties
%     x(8) = w_air      Air flow rate in air cooler
%     x(9) = Pf         Fan duty
%    x(10) = Tci_F      Water inlet temperature in furnace
%    x(11) = w_B        Total flow rate in hot water loop
%
%
% Helge Smedsrud
% IKP, NTNU
% 2008
% =========================================================================
%
clc
clear all
close all
%
% =========================================================================
%                                CONSTANTS
% =========================================================================
%
% Pressures [bara]
p1B = 15;                           % Brobekk side discharge pressure
p1V = 15;                           % Viken side discharge pressure
dp_f = 1000;                        % AC fan pressure difference (NB! [Pa])
%
% Densities [kg/m3]
rho_h = 912.892;                    % Hot water (p = 1500000 Pa, T = 155 °C)
rho_c = 959.049;                    % Cold water (p = 1500000 Pa, T = 100 °C)
rho_air_c = 1.161;                  % Cold air (p = 100000 Pa, T = 26.85 °C)
%
% Specific heat capacities [J/kg*K]
cph = 4321.56;                      % Hot water (p = 1500000 Pa, T = 155 °C)
cpc = 4213.84;                      % Cold water (p = 1500000 Pa, T = 100 °C)
cp_air = 1018.5;                    % Air (p = 100000 Pa, T = 126.85 °C)
%
% UA numbers [W/K]
UA_F = 52250.59846095464;           % Furnace
%
% Flow coefficients [m3/h*bar]
Kv_HB = 83.302787347224381;         % Brobekk side bypass valve
Kv_AC = 260.3212104600763;          % Brobekk side air cooler valve
Kv_HE = 260.3212104600763;          % Brobekk side heat exchanger valve
%
% Temperatures [°C]
Thi_F = 1000;                       % Furnace flue gas inlet
Tci_air = 10;                       % AC air inlet
%
% Flow rates [kg/s]
w_fg = 19.889359351988219;          % Furnace flue gas
%
% Areas [m2]
A_f = pi;                           % Fan duct
A_HE = 74.6;                        % Heat transfer area in heat exchanger
A_AC = A_HE;                        % Heat transfer area in air cooler
%
% Efficiencies [-]
z = 0.6;                            % Fraction of air velocity remaining in
                                    % diffusor area of AC fan
eta_p = 0.9;                        % Pump
eta_f = 0.9;                        % Fan
%
% Viken side inlet temperature [°C]
Tci_low = 65;                       % Minimum value
Tci_high = 90;                      % Maximum value
Tci_step = 0.1;                     % dT for the optimization 
%
% Viken side flow rate [kg/s]
w_V_low = 500000/3600;              % Minimum value
w_V_high = 900000/3600;             % Maximum value
w_V_step = (w_V_high-w_V_low)*Tci_step/...
           (Tci_high-Tci_low);      % dw for the optimization
%
% =========================================================================
%                            PROBLEM DEFINITION
% =========================================================================
%
% Initial values
%For Tci = 65.00 °C, and w_V = 500 t/h
x_0 = [0.925925925925397;0;1;0.279993940715972;1;0.067755925149853;
       0.001008737989590;0;0;1;1];
%
% Linear constraints (A needs to be non-empty in TOMLAB)
A = zeros(1,length(x_0));
A(1,1) = 1;
A = sparse(A);
%
% Lower and upper bounds on linear constraints
b_L = 0;
b_U = 1;
%
% Lower and upper bounds on nonlinear constraints
c_L = [90/135;90/149;110/126;0;0;0;0;0];
c_U = [1;1;1;0;0;0;0;0];
%
% Lower and upper bounds on x
x_L = [0;0;0.55;0;0.25;0;0;0;0;1;1];
x_U = [1;1;1;1;1;1;1;1;1;1;1];
%
% Setting up conAssign
Name = 'Brobekk incineration plant';
f = 'Brobekk_fun';
c = 'Brobekk_nonlincon';
g = [];
H = [];
HessPattern = [];
pSepFunc = [];
fLowBnd = 0;
dc = [];
d2c = [];
ConsPattern = [];
%
% Updating the TOMLAB problem structure
Prob = conAssign(f,g,H,HessPattern,x_L,x_U,Name,x_0,pSepFunc,...
                 fLowBnd,A,b_L,b_U,c,dc,d2c,ConsPattern,c_L,c_U);
%
% Cost vector
cost = [0;0;0;0;0;2;1;0;2;0;0];
%
% Scaling x and C using their largest allowed values
scx = [1;1;1;1;1;75000;75000;34.135734195475870;30000;126;250000/3600];
scC = [135;149;126];
%
% Adding the constants to Prob
Prob.p1B = p1B;             Prob.p1V = p1V;         Prob.rho_h = rho_h;
Prob.rho_c = rho_c;         Prob.Kv_HB = Kv_HB;     Prob.Kv_AC = Kv_AC;
Prob.Kv_HE = Kv_HE;         Prob.Tci_air = Tci_air; Prob.Thi_F = Thi_F;
Prob.cost = cost;           Prob.w_fg = w_fg;       Prob.rho_air_c = rho_air_c;
Prob.A_f = A_f;             Prob.z = z;             Prob.dp_f = dp_f;
Prob.eta_p = eta_p;         Prob.eta_f = eta_f;     Prob.scx = scx;
Prob.scC = scC;             Prob.UA_F = UA_F;       Prob.A_HE = A_HE;
Prob.A_AC = A_AC;           Prob.cp_air = cp_air;   Prob.cph = cph;
Prob.cpc = cpc;
%
% Creating vectors that will span the disturbance space
w_V_ref = w_V_low:w_V_step:w_V_high;
Tci_ref = Tci_low:Tci_step:Tci_high;
%
% Identifying the number of DOFs and disturbances
m = length(w_V_ref);
n = length(Tci_ref);
o = length(x_0);
%
% Preallocating matrices
xtflg = ones(m,n);
cst = zeros(m,n);
X = -1*ones(m,n,o);
%
% Additional user choices
PriorityDirection = 'T';
IterDetail = [];
RandFactor = 0.01;
%
% =========================================================================
%                               OPTIMIZATION
% =========================================================================
%
tic
% Go through the disturbance space spanned by w_V_ref and Tci_ref
for h = 1:m
    Prob.w_V = w_V_ref(h);
    disp(w_V_ref(h))
    for i = 1:n
        Prob.Tci = Tci_ref(i);
        %
        % Defining the scaling vector for Ceq
        Prob.scCeq = [16205850;250000/3600;w_V_ref(h);126;30000];
        %
        % Invoking Minos
        Result = tomRun('Minos',Prob,IterDetail);
        %
        % Convergence. The solution is used as x_0 for the next point
        if Result.ExitFlag == 0 && Result.Inform == 0
            xtflg(h,i) = Result.ExitFlag;
            cst(h,i) = cost'*(Result.x_k.*scx);
            X(h,i,:) = Result.x_k;
            Prob.x_0 = Result.x_k;
        %
        % No convergence. The routine will attempt to find a close point 
        % with convergence, randomly permute the {x} at this location a
        % small value, and then use this set as x_0 for the next point
        elseif Result.ExitFlag ~= 0 || Result.Inform ~= 0
            %
            % Problem in first column; search among the previous (1,i) 
            % points for feasible x_0 sets
            if h == 1 && i < n
                GG = xtflg(h,:);
                [uu,vv] = find(GG == 0);
                if isempty([uu,vv]) == 0
                    v = zeros(length(vv),1);
                    for k = 1:length(vv)
                        v(k) = abs(i-vv(k));
                    end
                    VV = find(min(v));
                    for j = 1:o
                        x0(j) = X(h,VV,j)+(1-X(h,VV,j))*rand*RandFactor;
                    end
                    Prob.x_0 = x0;
                else
                    disp('No zero exitflags in first column. Optimization terminated.')
                    return
                end
            %
            % Problem in second, third, fourth column etc.
            elseif h > 1 && i < n
                % Is it possible to use the set from (h-1,i+1)?
                if xtflg(h-1,i+1) == 0
                    x0 = X(h-1,i+1,:);
                    Prob.x_0 = x0;
                % Scanning the nearest disturbance space for feasible points
                elseif xtflg(h-1,i+1) ~= 0
                    if h < 4
                        GG = xtflg(1:h,:);
                    elseif h >= 4
                        GG = xtflg(h-2:h,:);
                    end
                    [uu,vv] = find(GG == 0);
                    lu = length(uu);
                    lv = length(vv);
                    if isempty([uu,vv]) == 0
                        u = zeros(lu,1);
                        v = zeros(lv,1);
                        for k = 1:lu
                            u(k) = abs(h-uu(k));
                            v(k) = abs(i-vv(k));
                        end
                        UU = find(min(u));
                        VV = find(min(v));
                        if UU < VV
                            for j = 1:o
                                x0(j) = X(UU,i,j)+(1-X(UU,i,j))*rand*RandFactor;
                            end
                            Prob.x_0 = x0;
                        elseif UU > VV
                            for j = 1:o
                                x0(j) = X(h,VV,j)+(1-X(h,VV,j))*rand*RandFactor;
                            end
                            Prob.x_0 = x0;
                        elseif UU == VV
                            if PriorityDirection == 'T';
                                for j = 1:o
                                    x0(j) = X(h,VV,j)+(1-X(h,VV,j))*rand*RandFactor;
                                end
                                Prob.x_0 = x0;
                            elseif PriorityDirection == 'w';
                                for j = 1:o
                                    x0(j) = X(UU,i,j)+(1-X(UU,i,j))*rand*RandFactor;
                                end
                                Prob.x_0 = x0;
                            end
                        end
                    end
                elseif isempty([uu,vv]) == 1
                    disp('No zero exitflags found. Optimization terminated.')
                    return
                end
            end
        end
        %
        % Check if the current point is at the end of the primary 
        % disturbance vector. If true, the routine wil "reset" x_0
        if i == n
            if h == 1 && n > 1
                if xtflg(h,1) == 0
                    x0 = X(h,1,:);
                    Prob.x_0 = x0;
                else
                    disp('Non-zero exitflag at (h,i).')
                    disp('Optimization terminated.')
                    return
                end
            elseif h > 1
                if xtflg(h,1) == 0 && xtflg(h-1,1) == 0
                    for j = 1:o
                        xdiff = X(h,1,j) - X(h-1,1,j);
                        x0(j) = X(h,1,j) + xdiff;
                    end
                    Prob.x_0 = x0;
                else
                    disp('Non-zero exitflags at (h,i) and (h-1,1).') 
                    disp('Optimization terminated.')
                    return
                end
            end
        end
    end
end
toc
%
% =========================================================================
%                                 PLOTTING
% =========================================================================
%
% If optimizing over a portion of the disturbance space larger than one
% point
if length(Tci_ref) > 1 && length(Tci_ref) == length(w_V_ref)
%  
% Generating the surface plots
for s = 1:o-2
    figure(s)
    zplt = X(:,:,s)*scx(s);
    surf(Tci_ref,w_V_ref,zplt)
    title(['Surface plot of optimization variable ' num2str(s)])
    xlabel('Cold water inlet temperature [°C]')
    ylabel('Cold water flow rate [kg s^-^1]')
    axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) ...
          x_L(s)*scx(s) x_U(s)*scx(s)])
end
%
% Adding z labels
figure(1)
zlabel('Brobekk side bypass valve position')
figure(2)
zlabel('Brobekk side air cooler valve position')
figure(3)
zlabel('Brobekk side heat exchanger valve position')
figure(4)
zlabel('Viken side bypass valve position')
figure(5)
zlabel('Viken side heat exchanger valve position')
figure(6)
zlabel('Brobekk side pump duty [W]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 0 15000])
figure(7)
zlabel('Viken side pump duty [W]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 0 1500])
figure(8)
zlabel('Flow rate of air in air cooler [kg s^-^1]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 0 20])
figure(9)
zlabel('Brobekk side air cooler fan duty [W]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 0 20000])
%
% Invoking sub-plot for calculation of additional variables
zzplt = -1*ones(length(w_V_ref),length(Tci_ref),4);
for i = 1:length(w_V_ref)
    for j = 1:length(Tci_ref)
        for k = 1:4
            zzplt(i,j,k) = Brobekk_calc(X(i,j,:),Prob,w_V_ref(i),Tci_ref(j),k);
        end
    end
end
%
% Additional plots
figure(o-1)
surf(Tci_ref,w_V_ref,zzplt(:,:,1))
title('Surface plot of T_V_i_k_e_n')
xlabel('Cold water inlet temperature [°C]')
ylabel('Cold water flow rate [kg s^-^1]')
zlabel('Temperature towards Viken/Oslo [°C]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 90 135])
%
figure(o)
surf(Tci_ref,w_V_ref,zzplt(:,:,2))
title('Surface plot of HE primary side outlet temperature')
xlabel('Cold water inlet temperature [°C]')
ylabel('Cold water flow rate [kg s^-^1]')
zlabel('Outlet temperature in HE, primary side [°C]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 110 126])
%
figure(o+1)
surf(Tci_ref,w_V_ref,zzplt(:,:,3))
title('Surface plot of HE secondary side outlet temperature')
xlabel('Cold water inlet temperature [°C]')
ylabel('Cold water flow rate [kg s^-^1]')
zlabel('Outlet temperature in HE, secondary side [°C]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 90 149])
%
figure(o+2)
surf(Tci_ref,w_V_ref,cst)
title('Surface plot of the operational cost')
xlabel('Cold water inlet temperature [°C]')
ylabel('Cold water flow rate [kg s^-^1]')
zlabel('Cost function value')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 0 50000])
%
figure(o+3)
surf(Tci_ref,w_V_ref,zzplt(:,:,4))
title('Surface plot of AC primary side outlet temperature')
xlabel('Cold water inlet temperature [°C]')
ylabel('Cold water flow rate [kg s^-^1]')
zlabel('Outlet temperature in AC, primary side [°C]')
axis([Tci_ref(1) Tci_ref(end) w_V_ref(1) w_V_ref(end) 110 180])
%
% Saving all parameters and variables in the entire script
save('All_data.mat')
%
% If optimizing on one single point in the disturbance space
elseif length(Tci_ref) == 1 && length(Tci_ref) == length(w_V_ref)
    format long
    Exitflag = Result.ExitFlag
    Cost = cst
    uHB = Result.x_k(1)*scx(1)
    uAC = Result.x_k(2)*scx(2)
    uHE = Result.x_k(3)*scx(3)
    uCB = Result.x_k(4)*scx(4)
    uCI = Result.x_k(5)*scx(5)
    PB = Result.x_k(6)*scx(6)
    PV = Result.x_k(7)*scx(7)
    w_air = Result.x_k(8)*scx(8)
    wB = Result.x_k(11)*scx(11)
    Pf = Result.x_k(9)*scx(9)
    Tci_F = Result.x_k(10)*scx(10)
    T_viken = Brobekk_calc(X,Prob,w_V_ref,Tci_ref,1)
    Tho_HE = Brobekk_calc(X,Prob,w_V_ref,Tci_ref,2)
    Tco_HE = Brobekk_calc(X,Prob,w_V_ref,Tci_ref,3)
    qF = Brobekk_calc(X,Prob,w_V_ref,Tci_ref,4)
    qHE = Brobekk_calc(X,Prob,w_V_ref,Tci_ref,5)
    qAC = Brobekk_calc(X,Prob,w_V_ref,Tci_ref,6)
    dq = qF - qHE - qAC
    Result.x_k
    format short
end