clear all
global Tk
Tk = 273.15;
open Ph.fig
flag = input('Design (1) or operation (2): ');
if flag==1
  d=data;
  Ths = Tk+25;
  Tcs = Tk-12;
  DTmin = 5;
  n = 1;
  [x, sc, d] = initialize(Ths, Tcs, DTmin, n, d);
  x_0 = x;
  Qc = 20e3;
elseif flag==2
  load data
  load optimal_design
  load scaling
  Ths = Tk+25;
  Tcs = Tk-12;
  x_0=x;
  con = x(1:23);
  vap = x(27:33);
  val = x(24:26);
  com = x(34:35);
  Acon = sum(con(6:10));
  Avap = vap(2);
  Qc = 20e3;  
end

% Condenser
Xcon_L = [zeros(5,1); zeros(5,1); -inf*ones(5,1); zeros(4,1); zeros(2,1); 0; 0]./sc.con.x;
Xcon_U = [inf*ones(5,1); inf*ones(5,1); inf*ones(5,1); inf*ones(4,1); inf*ones(2,1); 1; inf]./sc.con.x;

% Valve
Xval_L = [zeros(2,1); -inf]./sc.val.x;
Xval_U = [1; inf; inf]./sc.val.x;

% Vaporizer
Xvap_L = [0; 0; -inf; zeros(2,1); 0; 0]./sc.vap.x;
Xvap_U = [inf; inf; inf; inf*ones(2,1); inf; inf]./sc.vap.x;

% Compressor
Xcom_L = [0; 0]./sc.com.x;
Xcom_U = [inf; inf]./sc.com.x;

sc.x = [sc.con.x; sc.val.x; sc.vap.x; sc.com.x];
sc.res = [sc.con.res; sc.val.res; sc.vap.res; sc.com.res];

u = [Tcs; Ths];
aux.u = u;
aux.sc = sc;
aux.d = d;
% Use TOMLAB
 
res = nonlcon(x_0./sc.x, aux);
f = 'objfun';
g = [];
H = [];
HessPattern = [];
x_L = [Xcon_L; Xval_L; Xvap_L; Xcom_L];
x_U = [Xcon_U; Xval_U; Xvap_U; Xcom_U];
Name = 'Ammonia cycle';
pSepFunc = [];
fLowBnd = [];

if flag == 1
  A = xlsread('A.xls', 'design');
  b_L = [
         Qc;                    % Heat load in vaporizer
         (Ths+DTmin)*ones(5,1); % Delta T in condenser
         (Tcs-inf);             % Delta T in vaporizer
         ];
  b_U = [
         Qc;                    % Heat load in vaporizer
         (Ths+inf)*ones(5,1);   % Delta T in condenser
         (Tcs-DTmin);           % Delta T in vaporizer
         ];
elseif flag == 2
  A = xlsread('A.xls', 'operation');
  b_L = [
         Qc;                    % Heat load in vaporizer
         Acon;                  % Condenser area
         Avap;                  % Vaporizer area
         -inf;                  % Only used in operation
         -inf;                  % Only used in operation
         ];
  b_U = [
         Qc;                    % Heat load in vaporizer
         Acon;                  % Condenser area
         Avap;                  % Vaporizer area
         inf;                   % Only used in operation
         inf;                   % Only used in operation
         ];
end

SCALER = [];
for i=1:size(A,1)
  SCALER = [SCALER; sc.x'];
end
A=A.*SCALER;
c = 'nonlcon';

dc = [];
d2c = [];
ConsPattern = [];
c_L = zeros(size(res));
c_U = zeros(size(res));
x_min = [];
x_max = [];
f_opt = [];
x_opt = [];
Prob = ProbDef;
Prob = conAssign(f, g, H, HessPattern, x_L, x_U, Name, x_0./sc.x, pSepFunc, fLowBnd, A, b_L, b_U, c, dc, d2c, ConsPattern, c_L, c_U, x_min, x_max, f_opt, x_opt);
Prob.d = d;
Prob.sc = sc;

Prob.u = u;
Result = tomRun('minos', Prob);
disp(Result.ExitText)
x = Result.x_k.*sc.x;

if flag==1
  save optimal_design x
  save data d 
  save scaling sc
elseif flag==2
  save optimal_operation x
  save conditions u
end
xcon = x(1:23);
xvap = x(27:33);
xval = x(24:26);
xcom = x(34:35);

results
