Simultaneous synthesis of a heat exchanger network with multiple utilities using utility substages

Henryh · April 16, 2025, 11:47am

I try to model the “Simultaneous synthesis of a heat exchanger network with multiple utilities using utility substages
Jonggeol Na, Jaeheum Jung, Chansaem Park, Chonghun Han∗
School of Chemical and Biological Engineering, Seoul National University, Gwanak-ro 1, Gwanak-gu, Seoul 151-744, South Korea” paper of 28 equations motivating from the SYNHEAT posted code online in GAMS public, my code trying to model above paper MINLP is: * ==============================

SCALARS
==============================

SCALAR
NOK ‘number of stages in superstructure’ /3/
NOM ‘number of hot utilities’ /3/
NON ‘number of cold utilities’ /1/

EMAT ‘minimum approach temperature’ /10/
MT ‘total heat contents of stream i or j’ /1000/
theta ‘the maximum number of multiple utilities’ /2/
beta_cost ‘exponent for heat exchanger area cost’ /1/
betacu_cost ‘exponent for heat exchanger area cost for coldutility’ /1/
betahu_cost ‘exponent for heat exchanger area cost for hot utility’ /1/;
==============================
SET DEFINITIONS
==============================

$onText
HP = {i|i is a hot process stream}
CP = {j|j is a cold process stream}
HU = {m|m is a hot utility, m = 1,. . . NOM}
CU = {n|n is a cold utility, n = 1,. . . NON}
HU′ = {m|m is a hot utility, m = 1,. . . NOM + 1}
CU′ = {n|n is a cold utility, n = 1,. . . NON + 1}
ST = {k|k is a stage in the superstructure, k = 1,. . . NOK}
ST′ = {k|k is a stage in the superstructure, k = 1,. . . NOK + 1}
$offText

SET
i ‘hot process streams’ /12/ ,
j ‘cold process streams’ /11/ ,

USER NEED TO UPDATE HERE BASED ON SCALAR ABOVE of NOK, NOM, NON
m 'hot utilities - NOM +1 ’ /14/ ,
n ‘cold utilities - NON + 1’ /12/ ,
k ‘temperature locations - NOK +1’ /14/;

*USER NEED TO UPDATE HERE BASED ON SCALAR ABOVE of NOK, NOM, NON

NOMset ‘HU = {m|m is a hot utility, m = 1,. . . NOM}’ /1*2/,
NONset ‘CU = {n|n is a cold utility, n = 1,. . . NON}’ /1*2/,
NOMp1set ‘HU′ = {m|m is a hot utility, m = 1,. . . NOM + 1}’ /1*3/,
NONp1set ‘CU′ = {n|n is a cold utility, n = 1,. . . NON + 1}’ /1*3/,
NOKset ‘ST = {k|k is a stage in the superstructure, k = 1,. . . NOK}’ /1*2/,
NOKp1set ‘ST′ = {k|k is a stage in the superstructure, k = 1,. . . NOK + 1}’ /1*3/;

Set
HP(i)
CP(j)
HU(m)
CU(n)
ST(k)
HUP(m)
CUP(n)
STP(k);

STP(k) = yes$(ord(k) <= NOK+1);
HUP(m) = yes$(ord(m) <= NOM+1);
CUP(n) = yes$(ord(n) <= NON+1);

Singleton Set
firstk(k) ‘first k temperature location’
lastk(k) ‘last k temperature location’

firstm(m) ‘first m temperature location’
lastm(m) ‘last m temperature location’

firstn(n) 'first n temperature location'

lastn(n) ‘last n temperature location’;

ST(k) = yes$(ord(k) < card(k));
firstk(k) = yes$(ord(k) = 1);
lastk(k) = yes$(ord(k) = card(k));

HU(m) = yes$(ord(m) < card(m));
firstm(m) = yes$(ord(m) = 1);
lastm(m) = yes$(ord(m) = card(m));

CU(n) = yes$(ord(n) < card(n));
firstn(n) = yes$(ord(n) = 1);
lastn(n) = yes$(ord(n) = card(n));

==============================
PARAMETER DEFINITIONS
==============================

PARAMETER
A(i,j,k) ‘heat exchanger area between hot stream i and cold stream j at stage k’
Acu(i,k,n) ‘heat exchanger area of cold utility n and hot stream i at stage k’
Ahu(j,k,m) ‘heat exchanger area of hot utility m and cold stream j at stage k’
C(i,j) ‘area cost coefficient for heat exchanger between hot stream i and cold stream j’
CC(i,n) ‘area cost coefficient for cold utility n’
CCU(n) ‘per unit cost for cold utility n’
CF(i,j) ‘fixed cost for process stream heat exchanger between hot stream i and cold stream j’
CFC(i,n) ‘fixed cost for cold utility heat exchanger for cold utility n’
CFH(j,m) ‘fixed cost for hot utility heat exchanger for hot utility m’
CH(j,m) ‘area cost coefficient for hot utility heat exchangers’
CHU(m) ‘per unit cost for hot utility’

EMAT
Fhp(i)        'heat capacity flow rate of hot stream i'
Fcp(j)        'heat capacity flow rate of cold stream j'
hhp(i)        'heat transfer coefficient for hot stream i'
hcp(j)        'heat transfer coefficient for cold stream j'
hhu(m)        'heat transfer coefficient for hot utility m'
hcu(n)        'heat transfer coefficient for cold utility n'

LMTD(i,j,k) 'Chen’s log-mean temperature difference for hot stream i and cold stream j at stage k'
LMTDP(i,j,k) 'Paterson’s log-mean temperature difference for hot stream i and cold stream j at stage k'


TINcu(n)    'inlet temperature of cold utility'
TINhu(m)    'inlet temperature of hot utility'
TINhp(i)      'inlet temperature of hot stream i'
TINcp(j)      'inlet temperature of cold stream j'
TOUTcu(n)   'outlet temperature of cold utility'
TOUThu(m)   'outlet temperature of hot utility'
TOUThp(i)     'outlet temperature of hot stream i' 
TOUTcp(j)     'outlet temperature of cold stream j' ;

*PARAMETER

gammaij(i,j) ‘gamma for i-j matches’
gammajm(j,m) ‘gamma for j-m utility matches’
gammain(i,n) ‘gamma for i-n utility matches’;
==============================
BINARY VARIABLES
==============================

BINARY VARIABLE
z(i,j,k) ‘binary variable denoting existence of heat exchanger between hot stream i and cold stream j at stage k’
zcu(i,k,n) ‘binary variable denoting existence of heat exchanger between hot stream i and cold utility n at stage k’
zhu(j,k,m) ‘binary variable denoting existence of heat exchanger between hot utility m and cold stream j at stage k’ ;

==============================
POSITIVE VARIABLES
==============================

Positive Variable
q(i,j,k) ‘heat exchanged between process stream i and j at stage k’
qcu(i,k,n) ‘heat exchanged between hot utility m and cold stream j at stage k’
qhu(j,k,m) ‘heat exchanged between hot stream i and cold utility n at stage k’

tc(j,k,m)   'temperature of cold stream j at stage k and utility substage m'
th(i,k,n)   'temperature of hot stream i at stage k and utility substage n'

dtcu(i,k,n) 'temperature approach for matching hot stream i and cold utility n at stage k'
dthu(j,k,m) 'temperature approach for matching hot utility m and cold stream j at stage k'
dtl(i,j,k)  'left: temperature approach for matching stream i and j at stage k'
dtr(i,j,k)  'right: temperature approach for matching stream i and j at stage k';

Variable cost ‘hen and utility cost’;

Equation
eq1 ‘Overall heat balance for hot process stream i in stage k’
eq2 ‘Overall heat balance for cold process stream j in stage k’

eq3 ''
eq4 ''

eq5 ''
eq6 ''

eq7 ''
eq8 ''
eq9 ''
eq10 ''

eq11 ''
eq12 ''
eq13 ''
eq14 ''

eq15 ''
eq16 ''
eq17 ''

eq18m 'optionally'
eq18n 'optionally'

eq19
eq20
eq21
eq22
eq23
eq24

eq25a
eq25b
eq25c
eq25d


obj
;

eq1(i).. (TINhp(i) - TOUThp(i))*Fhp(i) =e= sum((j,ST), q(i,j,ST)) + sum((STP,CU),qcu(i,STP,CU));
eq2(j).. (TINcp(j) - TOUTcp(j))*Fcp(j) =e= sum((i,ST), q(i,j,ST)) + sum((STP,HU),qhu(j,STP,HU));

eq3(i,k)$ST(k).. (th(i,k,lastn) - th(i,k+1,firstn))*Fhp(i) =e= sum(j, q(i,j,k));
eq4(j,k)$ST(k).. (tc(j,k,lastm) - tc(j,k+1,firstm))*Fcp(j) =e= sum(i, q(i,j,k));

eq5(i,k,n)$ST(k).. (th(i,k,n) - th(i,k,n+1))*Fhp(i) =e= qcu(i,k,n);
eq6(j,k,m)$ST(k).. (tc(j,k,m) - tc(j,k,m+1))*Fcp(j) =e= qhu(j,k,m);

eq7(i).. th(i,firstk,firstn) =e= TINhp(i);
eq8(j).. tc(j,firstk,firstm) =e= TOUTcp(j);
eq9(i).. th(i,lastk,lastn) =e= TOUThp(i);
eq10(j).. tc(j,lastk,lastm) =e= TINcp(j);

eq11(i,k,n)$ST(k).. th(i,k,n) =g= th(i,k,n+1);
eq12(j,k,m)$ST(k).. tc(j,k,m) =g= tc(j,k,m+1);
eq13(i,k)$ST(k).. th(i,k,lastn) =g= th(i,k+1,firstn);
eq14(j,k)$ST(k).. tc(j,k,lastm) =g= tc(j,k+1,firstm);

ech(i) = Fhp(i)(TINhp(i) - TOUThp(i));
ecc(j) = Fcp(j)(TOUTcp(j) - TOUTcp(j));

eq15(i,j,k)$ST(k).. q(i,j,k) - min(Fhp(i)(TINhp(i) - TOUThp(i)), Fcp(j)(TOUTcp(j) - TOUTcp(j)))z(i,j,k) =l= 0;
eq16(i,k,n)$ST(k).. qcu(i,k,n) - (Fhp(i)(TINhp(i) - TOUThp(i)))zcu(i,k,n) =l= 0;
eq17(j,k,m)$ST(k).. qhu(j,k,m) - (Fcp(j)(TOUTcp(j) - TOUTcp(j)))*zhu(j,k,m) =l= 0;

eq18m(j,k)$ST(k).. sum(m,zhu(j,k,m)) =l= theta;
eq18n(i,k)$ST(k).. sum(n,zcu(i,k,n)) =l= theta;

*gammaij(i,j) = max(0,TINcp(j) - TINhp(i), TINcp(j) - TOUThp(i),TOUTcp(j) - TINhp(i), TOUTcp(j) - TOUThp(i));
*gammajm(j,m) = max(0,TINcu(m) - TINcp(j), TINcu(m) - TOUTcp(j),TOUTcu(m) - TINcp(j), TOUTcu(m) - TOUTcp(j));
*gammain(i,n) = max(0,TINhu(n) - TINhp(i), TINhu(n) - TOUThp(i),TOUThu(n) - TINhp(i), TOUThu(n) - TOUThp(i));

eq19(i,j,k)$ST(k)..dtl(i,j,k) =l= th(i,k,lastn) - tc(j,k,lastm) + max(0,TINcp(j) - TINhp(i), TINcp(j) - TOUThp(i),TOUTcp(j) - TINhp(i), TOUTcp(j) - TOUThp(i))(1 - z(i,j,k));
eq20(i,j,k)$ST(k)..dtr(i,j,k) =l= th(i,k+1,firstn) - tc(j,k+1,firstm) + max(0,TINcp(j) - TINhp(i), TINcp(j) - TOUThp(i),TOUTcp(j) - TINhp(i), TOUTcp(j) - TOUThp(i))(1 - z(i,j,k));

eq21(j,k,m,n)$ST(k)..dthu(j,k,m) =l= TINhu(m) - tc(j,k,m) + max(0,TINcu(n) - TINcp(j), TINcu(n) - TOUTcp(j),TOUTcu(n) - TINcp(j), TOUTcu(n) - TOUTcp(j))(1-zhu(j,k,m));
eq22(j,k,m,n)$ST(k)..dthu(j,k,m+1) =l= TOUThu(m) - tc(j,k,m+1) + max(0,TINcu(n) - TINcp(j), TINcu(n) - TOUTcp(j),TOUTcu(n) - TINcp(j), TOUTcu(n) - TOUTcp(j))(1-zhu(j,k,m));

eq23(i,k,n,m)$ST(k)..dtcu(i,k,n) =l= th(i,k,n) - TOUTcu(n) + max(0,TINhu(m) - TINhp(i), TINhu(m) - TOUThp(i),TOUThu(m) - TINhp(i), TOUThu(m) - TOUThp(i))(1-zcu(i,k,n));
eq24(i,k,n,m)$ST(k)..dtcu(i,k,n+1) =l= th(i,k,n+1) - TINcu(n) + max(0,TINhu(m) - TINhp(i), TINhu(m) - TOUThp(i),TOUThu(m) - TINhp(i), TOUThu(m) - TOUThp(i))(1-zcu(i,k,n));

eq25a(i,j,k)..dtl(i,j,k) =l= EMAT;
eq25b(i,j,k)..dtr(i,j,k) =l= EMAT;
eq25c(j,k,m)..dthu(j,k,m) =l= EMAT;
eq25d(i,k,n)..dtcu(i,k,n) =l= EMAT;

obj..cost =e=

(
sum((i,k,n),CCU(n)*qcu(i,k,n))

+ sum((j,k,m),CHU(m)*qhu(j,k,m))

)

+(
sum(
(i,j,k), CF(i,j) * z(i,j,k))
+ sum((i,k,n), CFC(i,n) * zcu(i,k,n))
+ sum((j,k,m),CFH(j,m) * zhu(j,k,m))

)

          sum((i,j,k),
          C(i,j) * (
              (
              (q(i,j,k) *((1/hhp(i)) + (1/hcp(j)))
              /
              ((dtl(i,j,k)*dtr(i,j,k+1) * (dtl(i,j,k) + dtr(i,j,k+1))/2
                  + 1e-6)**0.33333) + 1e-6) + 1e-6
              )**beta_cost
              )
          )

sum((i,k,n),
CC(i,n)* (
(qcu(i,k,n)*((1/hhp(i)) + (1/hcu(n)))
/ (((TOUThp(i)-TINcu(n))dtcu(i,k,n)(TOUThp(i) - TINcu(n) + dtcu(i,k,n))/2
+ 1e-6)**0.33333) + 1e-6)**betacu_cost
)
)
sum((j,k,m),
CH(j,m)* (
(qhu(j,k,m)*((1/hcp(j)) + 1/hhu(m)))
/ (((TINhu(m) - TOUTcp(j))dthu(j,k,m)((TINhu(m) - TOUTcp(j) + dthu(j,k,m))/2)
+ 1e-6)**0.33333) + 1e-6)**betahu_cost
)

;

Process streams
hot
TINhp(‘1’) = 105; TOUThp(‘1’) = 25; Fhp(‘1’) = 10; hhp(‘1’) = 0.5;
TINhp(‘2’) = 185; TOUThp(‘2’) = 35; Fhp(‘2’) = 5; hhp(‘2’) = 0.5;
cold
TINcp(‘1’) = 25; TOUTcp(‘1’) = 185; Fcp(‘1’) = 7.5; hcp(‘1’) = 0.5;
*TINcp(‘2’) = 350; TOUTcp(‘2’) = 500; Fcp(‘2’) = 13; hcp(‘2’) = 1;
costs and coefficients
CHU(‘1’) = 160; TINhu(‘1’) = 210; TOUThu(‘1’) = 209; hhu(‘1’) = 5;
CHU(‘2’) = 110; TINhu(‘2’) = 160; TOUThu(‘3’) = 159; hhu(‘1’) = 5;
CHU(‘3’) = 50; TINhu(‘3’) = 130; TOUThu(‘3’) = 129; hhu(‘1’) = 5;

CCU(‘1’) = 10; TINcu(‘1’) = 5; TOUTcu(‘1’) = 6; hcu(‘1’) = 2.6;

CF(i,j) = 1;
C(i,j) = 150;
*aexp = 1;
*EMAT = 10;
EMAT = 20;

bounds

dtl.lo(i,j,k) = EMAT;
dtr.lo(i,j,k) = EMAT;
dthu.lo(j,k,m) = EMAT;
dtcu.lo(i,k,n) = EMAT;
th.up(i,k,n) = TINhp(i);
th.lo(i,k,n) = TOUThp(i);
tc.up(j,k,m) = TOUTcp(j);
tc.lo(j,k,m) = TINcp(j);

initialization
th.l(i,k,n) = TINhp(i);
tc.l(j,k,m) = TINcp(j);
dthu.l(j,k,m) = TOUThu(m) - TINcp(j);
dtcu.l(i,k,n) = TINhp(i) - TOUTcu(n);

ech(i) = Fhp(i)(TINhp(i) - TOUThp(i));
ecc(j) = Fcp(j)(TOUTcp(j) - TOUTcp(j));

*gammaij(i,j) = max(0,TINcp(j) - TINhp(i), TINcp(j) - TOUThp(i),TOUTcp(j) - TINhp(i), TOUTcp(j) - TOUThp(i));
*gammajm(j,m) = max(0,TINcu(m) - TINcp(j), TINcu(m) - TOUTcp(j),TOUTcu(m) - TINcp(j), TOUTcu(m) - TOUTcp(j));
*gammain(i,n) = max(0,TINhu(n) - TINhp(i), TINhu(n) - TOUThp(i),TOUThu(n) - TINhp(i), TOUThu(n) - TOUThp(i));

dtl.l(i,j,k) = TINhp(i) - TINcp(j);
dtr.l(i,j,k) = TINhp(i) - TINcp(j);
q.l(i,j,k)$ST(k) = min(Fhp(i)(TINhp(i) - TOUThp(i)),Fcp(j)(TOUTcp(j) - TOUTcp(j)));

Initialization for cost parameters
CFC(i,n) = 1;
CFH(j,m) = 1;
CC(i,n) = 1;
CH(j,m) = 1;

Model super/ all /;

option optCr = 0, limRow = 0, limCol = 0, solPrint = off, sysOut = off;

solve super using minlp minimizing cost;

================================ where I can not find any way to resolve problem of division by 0, error is described as: “— Starting compilation
— NEW_3_TEST.gms(461) 3 Mb
— Starting execution: elapsed 0:00:00.005[LST:465]
— NEW_3_TEST.gms(364) 4 Mb
— Generating MINLP model super[LST:465]
— NEW_3_TEST.gms(255) 6 Mb
*** Error at line 255: division by zero (0)[LST:469]
— NEW_3_TEST.gms(255) 6 Mb 1
*** Error at line 255: A constant in a nonlinear expression in equation obj evaluated to UNDF[LST:470]
— NEW_3_TEST.gms(255) 6 Mb 2
*** Error at line 255: rPower: FUNC DOMAIN: x**y, x < 0[LST:471]
— NEW_3_TEST.gms(366) 4 Mb 3
*** SOLVE aborted
— NEW_3_TEST.gms(366) 4 Mb 3
*** Status: Execution error(s)[LST:503]
— Job NEW_3_TEST.gms Stop” , the working SYNHEAT model is: ```
$title Simultaneous Optimization for HEN Synthesis (SYNHEAT,SEQ=117)

$onText
This model designs a heat exchanger network which operates at
minimal annual cost and satisfies heating and cooling require-
ments. The superstructure consists of two stages with eight
possible exchangers.

Yee, T F, and Grossmann, I E, Simultaneous Optimization of Models for
Heat Integration - Heat Exchanger Network Synthesis. Computers and
Chemical Engineering 14, 10 (1990), 1151-1184.

Keywords: mixed integer nonlinear programming, chemical engineering, heat exchanger
network
$offText

Set
i ‘hot streams’ / 12 /
j ‘cold streams’ / 12 /;

Scalar nok ‘number of stages in superstructure’ / 2 /;

Set
k ‘temperature locations nok + 1’ / 1*3 /
st(k) ‘stages’;

Singleton Set
firstK(k) ‘first temperature location’
lastK(k) ‘last temperature location’;

st(k) = yes$(ord(k) < card(k));
firstK(k) = yes$(ord(k) = 1);
lastK(k) = yes$(ord(k) = card(k));

Parameter
fh(i) ‘heat capacity flowrate of hot stream’
fc(j) ‘heat capacity flowrate of cold stream’
thin(i) ‘supply temp. of hot stream’
thout(i) ‘target temp. of hot stream’
tcin(j) ‘supply temp. of cold stream’
tcout(j) ‘target temp. of cold stream’
ech(i) ‘heat content hot I’
ecc(j) ‘heat content cold j’
hh(i) ‘stream-individual film coefficient hot I’
hc(j) ‘stream-individual film coefficient cold j’
hucost ‘cost of heating utility’
cucost ‘cost of cooling utility’
unitc ‘fixed charge for exchanger’
acoeff ‘area cost coefficient for exchangers’
hucoeff ‘area cost coefficient for heaters’
cucoeff ‘area cost coefficient for coolers’
aexp ‘cost exponent for exchangers’
hhu ‘stream-individual film coefficient hot utility’
hcu ‘stream-individual film coefficient cold utility’
thuin ‘inlet temperature hot utility’
thuout ‘outlet temperature hot utility’
tcuin ‘inlet temperature cold utility’
tcuout ‘outlet temperature cold utility’
gamma(i,j) ‘upper bound of driving force’
a(i,j,k) ‘area for exchanger for match ij in interval k (chen approx.)’
al(i,j,k) ‘area calculated with log mean’
acu(i) ‘area coolers’
ahu(j) ‘area heaters’
tmapp ‘minimum approach temperature’
costheat ‘cost of heating’
costcool ‘cost of cooling’
invcost ‘investment cost’;

Binary Variable
z(i,j,k)
zcu(i)
zhu(j);

Positive Variable
th(i,k) ‘temperature of hot stream i as it enters stage k’
tc(j,k) ‘temperature of cold stream j as it leaves stage k’
q(i,j,k) ‘energy exchanged between i and j in stage k’
qc(i) ‘energy exchanged between i and the cold utility’
qh(j) ‘energy exchanged between j and the hot utility’
dt(i,j,k) ‘approach between i and j at location k’
dtcu(i) ‘approach between i and the cold utility’
dthu(j) ‘approach between j and the hot utility’;

Variable cost ‘hen and utility cost’;

Equation
eh(i,k) ‘energy exchanged by hot stream i in stage k’
eqc(i) ‘energy exchanged by hot stream i with the cold utility’
teh(i) ‘total energy exchanged by hot stream I’
ec(j,k) ‘energy exchanged by cold stream j in stage k’
eqh(j) ‘energy exchanged by cold stream j with the hot utility’
tec(j) ‘total energy exchanged by cold stream j’
month(i,k) ‘monotonicity of th’
montc(j,k) ‘monotonicity of tc’
monthl(i) ‘monotonicity of th k = last’
montcf(j) ‘monotonicity of tc for k = 1’
tinh(i) ‘supply temperature of hot streams’
tinc(j) ‘supply temperature of cold streams’
logq(i,j,k) ‘logical constraints on q’
logqh(j) ‘logical constraints on qh(j)’
logqc(i) ‘logical constraints on qc(i)’
logdth(i,j,k) ‘logical constraints on dt at the hot end’
logdtc(i,j,k) ‘logical constraints on dt at the cold end’
logdtcu(i) ‘logical constraints on dtcu’
logdthu(j) ‘logical constraints on dthu’
obj ‘objective function’;

teh(i).. (thin(i) - thout(i))*fh(i) =e= sum((j,st), q(i,j,st)) + qc(i);
tec(j).. (tcout(j) - tcin(j))*fc(j) =e= sum((i,st), q(i,j,st)) + qh(j);

eh(i,k)$st(k).. fh(i)(th(i,k) - th(i,k+1)) =e= sum(j, q(i,j,k));
ec(j,k)$st(k).. fc(j)(tc(j,k) - tc(j,k+1)) =e= sum(i,q(i,j,k));

eqc(i).. fh(i)(th(i,lastK) - thout(i)) =e= qc(i);
eqh(j).. fc(j)(tcout(j) - tc(j,firstK)) =e= qh(j);

tinh(i).. thin(i) =e= th(i,firstK);
tinc(j).. tcin(j) =e= tc(j,lastK);

month(i,k)$st(k).. th(i,k) =g= th(i,k+1);
montc(j,k)$st(k).. tc(j,k) =g= tc(j,k+1);

monthl(i).. th(i,lastK) =g= thout(i);
montcf(j).. tcout(j) =g= tc(j,firstK);

logq(i,j,k)$st(k).. q(i,j,k) - min(ech(i), ecc(j))*z(i,j,k) =l= 0;

logqc(i).. qc(i) - ech(i)*zcu(i) =l= 0;
logqh(j).. qh(j) - ecc(j)*zhu(j) =l= 0;

logdth(i,j,k)$st(k)..
dt(i,j,k) =l= th(i,k) - tc(j,k) + gamma(i,j)*(1 - z(i,j,k));

logdtc(i,j,k)$st(k)..
dt(i,j,k+1) =l= th(i,k+1) - tc(j,k+1) + gamma(i,j)*(1 - z(i,j,k));

logdthu(j).. dthu(j) =l= (thuout - tc(j,firstK));
logdtcu(i).. dtcu(i) =l= th(i,lastK) - tcuout;

obj..cost =e= unitc*(sum((i,j,st),z(i,j,st))
+ sum(i,zcu(i)) + sum(j,zhu(j)))
+ acoeffsum((i,j,k),(q(i,j,k)((1/hh(i)) + (1/hc(j)))
/ (((dt(i,j,k)dt(i,j,k+1)(dt(i,j,k) + dt(i,j,k+1))/2
+ 1e-6)**0.33333) + 1e-6) + 1e-6)*aexp)
+ hucoeff(sum(j,(qh(j)((1/hc(j)) + 1/hhu))
/ (((thuin - tcout(j))dthu(j)((thuin - tcout(j) + dthu(j))/2)
+ 1e-6)**0.33333) + 1e-6)**aexp)
+ cucoeffsum(i,(qc(i)*((1/hh(i)) + (1/hcu))
/ (((thout(i)-tcuin)dtcu(i)(thout(i) - tcuin + dtcu(i))/2
+ 1e-6)**0.33333) + 1e-6)**aexp)
+ sum(j,qh(j)*hucost) + sum(i,qc(i)*cucost);

Process streams
hot
thin(‘1’) = 650; thout(‘1’) = 370; fh(‘1’) = 10; hh(‘1’) = 1;
thin(‘2’) = 590; thout(‘2’) = 370; fh(‘2’) = 20; hh(‘2’) = 1;
cold
tcin(‘1’) = 410; tcout(‘1’) = 650; fc(‘1’) = 15; hc(‘1’) = 1;
tcin(‘2’) = 350; tcout(‘2’) = 500; fc(‘2’) = 13; hc(‘2’) = 1;
costs and coefficients
hucost = 80; hucoeff = 150; thuin = 680; thuout = 680; hhu = 5;
cucost = 15; cucoeff = 150; tcuin = 300; tcuout = 320; hcu = 1;

unitc = 5500; acoeff = 150;
aexp = 1; tmapp = 10;

bounds
dt.lo(i,j,k) = tmapp;
dthu.lo(j) = tmapp;
dtcu.lo(i) = tmapp;
th.up(i,k) = thin(i);
th.lo(i,k) = thout(i);
tc.up(j,k) = tcout(j);
tc.lo(j,k) = tcin(j);
initialization
th.l(i,k) = thin(i);
tc.l(j,k) = tcin(j);
dthu.l(j) = thuout - tcin(j);
dtcu.l(i) = thin(i) - tcuout;
ech(i) = fh(i)(thin(i) - thout(i));
ecc(j) = fc(j)(tcout(j) - tcin(j));
gamma(i,j) = max(0,tcin(j) - thin(i), tcin(j) - thout(i),tcout(j) - thin(i), tcout(j) - thout(i));

dt.l(i,j,k) = thin(i) - tcin(j);
q.l(i,j,k)$st(k) = min(ech(i),ecc(j));

Model super/ all /;

option optCr = 0, limRow = 0, limCol = 0, solPrint = off, sysOut = off;

solve super using minlp minimizing cost;

Areas by chen approximation
a(i,j,k)$st(k) = q.l(i,j,k)((1/hh(i)) + (1/hc(j)))
/ (2/3sqrt(dt.l(i,j,k)dt.l(i,j,k+1))
+ 1/6(1e-8 + dt.l(i,j,k) + dt.l(i,j,k+1)));
Areas by log mean temperature
al(i,j,k)$st(k) = (q.l(i,j,k)*((1/hh(i)) + (1/hc(j))))
/ (dt.l(i,j,k)*dt.l(i,j,k+1)
* (dt.l(i,j,k) + dt.l(i,j,k+1))/2)**0.33333;

display a, al;

Areas of heaters and coolers
ahu(j) = (qh.l(j)*((1/hc(j)) + (1/hhu))/(((thuin-tcout(j))*dthu.l(j)
* ((thuin-tcout(j)+dthu.l(j)))/2) + 1e-6)**0.33333);

acu(i) = (qc.l(i)*((1/hh(i)) + (1/hcu))/(((thout(i) - tcuin)*dtcu.l(i)
* (thout(i) - tcuin + dtcu.l(i))/2 + 1e-6)**0.33333));

display acu, ahu;

Utility costs
costheat = sum(j,qh.l(j)*hucost);
costcool = sum(i,qc.l(i)*cucost);

display costheat, costcool;

Investment cost
invcost = cost.l - costheat - costcool;

display invcost;