This Mathematica script determines the optimal interstage pressures for a 4-stage propane refrigeration compression cycle (unequal cooling stages) that minimize total compressor work.
Thermodynamic property models: ideal-gas heat capacities (Cp) as well as the Soave–Redlich–Kwong (SRK) equation of state with departure functions (Hres, Sres) are used to compute real-gas enthalpy and entropy.
Stage-by-stage compression model: for each stage, I compute the isentropic outlet temperature (FindRoot on entropy balance) and the corresponding work input, accounting for both real-gas (SRK) and ideal-gas assumptions.
Objective function (WtotalSRK): total compressor work summed across all 4 stages, given a set of trial interstage pressures.
Custom Particle Swarm Optimization (PSO): an implementation (psoOptimize) with velocity/position updates, global/personal best tracking, and convergence history logging, is used to search for the pressure set (P2, P3, P4), which minimizes total work.
Two optimization runs: one using real-gas (SRK) properties, and one using ideal-gas properties — to quantify the impact of real-gas non-idealities on optimal staging and compressor power.
Convergence diagnostics: log-scale plots of best-loss-vs-iteration for both cases, enabling comparison of PSO convergence behavior and final optimized work/pressure results for either the real-gas case or the ideal-gas models.
Optimal Interstage Pressures for a 4-Stage Compressor with unequal Interstage Cooling -- Real Gas (SRK EoS) or IG, propane. Uses a custom Particle Swarm Optimization routine (psoOptimize)
ClearAll["Global`*"];
(* Propane physical constants *)
Rgas=8.314;(*J/(molK)*)Tc=369.8;(*K*)Pc=42.48*^5;(*Pa*)omega=0.152;Mmass=0.04410;(*kg/mol,molarmassofpropane*)
SRK parameters
aSRK=0.42748*Rgas^2*Tc^2/Pc;bSRK=0.08664*Rgas*Tc/Pc;mSRK=0.480+1.574*omega-0.176*omega^2;
In[]:=
alpha[T_]:=(1+mSRK*(1-Sqrt[T/Tc]))^2;alphaPrime[T_]:=-mSRK*(1+mSRK*(1-Sqrt[T/Tc]))/Sqrt[T*Tc];
Check
FullSimplify[D[alpha[T],T]-alphaPrime[T],Tc>0]
Out[]=
2.58061×-
-19
10
1.19015×
-17
10
T
Ideal-gas Cp0 correlation for propane: Cp0/R = A + B T + C T^2 (Smith, Van Ness, Abbott coefficients)
Acp=1.213;Bcp=28.785*^-3;Ccp=-8.824*^-6;Tref=300.0;(*K,referencestate*)Pref=1.0*^5;(*Pa(1bar)*)
Cp0overR[T_]:=Acp+Bcp*T+Ccp*T^2;
Ideal-gas enthalpy relative to Tref (J/mol)
Hig[T_?NumericQ]:=NIntegrate[Rgas*Cp0overR[Tt],{Tt,Tref,T}];
Ideal-gas entropy relative to (Tref, Pref) (J/mol/K)
Sig[T_?NumericQ,P_?NumericQ]:=NIntegrate[Rgas*Cp0overR[Tt]/Tt,{Tt,Tref,T}]-Rgas*Log[P/Pref];
SRK compressibility factor (vapor root - highest root - close to 0.9)
ZSRK[T_?NumericQ,P_?NumericQ]:=Module[{A,B,coeffs},A=aSRK*alpha[T]*P/(Rgas*T)^2;B=bSRK*P/(Rgas*T);(*Z^3-Z^2+(A-B-B^2)Z-AB=0*)coeffs={1.0,-1.0,(A-B-B^2),-A*B};z/.FindRoot[coeffs[[1]]*z^3+coeffs[[2]]*z^2+coeffs[[3]]*z+coeffs[[4]]==0,{z,0.9}]];
SRK residual (departure) enthalpy and entropy
Hres[T_?NumericQ,P_?NumericQ]:=Module[{Z,B,term},Z=ZSRK[T,P];B=bSRK*P/(Rgas*T);term=(aSRK*(alpha[T]-T*alphaPrime[T])/bSRK)*Log[Z/(Z+B)];Rgas*T*(Z-1)+term];Sres[T_?NumericQ,P_?NumericQ]:=Module[{Z,B},Z=ZSRK[T,P];B=bSRK*P/(Rgas*T);Rgas*Log[Z-B]-(aSRK*alphaPrime[T]/bSRK)*Log[Z/(Z+B)]];Hreal[T_?NumericQ,P_?NumericQ]:=Hig[T]+flagHres[T,P];Sreal[T_?NumericQ,P_?NumericQ]:=Sig[T,P]+flagSres[T,P];
Per-stage isentropic work (real gas, SRK). Given inlet (Tin,Pin), find Tout at Pout such that entropy is conserved (isentropic step), then return actual work = dH (J/mol) and Tout)
stageWork[TinStage_?NumericQ,PinStage_?NumericQ,PoutStage_?NumericQ]:=Module[{Sin,entropyResidual,Tguess,Tlo,Thi,flo,fhi,tries,Tout,wJmol},Sin=Sreal[TinStage,PinStage];entropyResidual[T_?NumericQ]:=Sreal[T,PoutStage]-Sin;(*FindTout*)Tout=T/.FindRoot[entropyResidual[T]==0,{T,331,330},MaxIterations->1000];wJmol=Hreal[Tout,PoutStage]-Hreal[TinStage,PinStage];{wJmol,Tout}];
Problem setup: 4 stages, unequal interstage cooling
P1=1.0*^5;(*Pa,suctionpressure(1bar)*)Pn=9.0*^5;(*Pa,finaldischargepressure(8bar)*)T1=300.0;(*K*)dT={5.0,8.0,12.0};Tin=Prepend[T1+Accumulate[dT],T1];(*{300,305,313,325}K*)Print["Stage suction temperatures (K): ",Tin];
Stage suction temperatures (K): {300.,305.,313.,325.}
Total work objective
WtotalSRK[{P2_?NumericQ,P3_?NumericQ,P4_?NumericQ}]:=Module[{Pstages,totalJPerMol,w},Pstages={P1,P2,P3,P4,Pn};totalJPerMol=0.0;Do[w=stageWork[Tin[[i]],Pstages[[i]],Pstages[[i+1]]][[1]];totalJPerMol+=w,{i,1,4}];(totalJPerMol/Mmass)(*J/kg*)];
Custom Particle Swarm Optimizer (same as ideal-gas notebooks)
clipMatrix[mat_,lbv_,ubv_]:=Transpose[MapThread[Clip[#1,{#2,#3}]&,{Transpose[mat],lbv,ubv}]];clampMatrix[mat_,vmaxv_]:=Transpose[MapThread[Clip[#1,{-#2,#2}]&,{Transpose[mat],vmaxv}]];
psoOptimize[lossFunc_,dim_,lb_,ub_,nParticles_:30,nIter_:100,wMax_:0.9,wMin_:0.4,c1_:1.42,c2_:1.42,vClampFrac_:0.2]:=Module[{particles,velocities,pbest,pbestScore,gbest,gbestScore,convergence={},r1,r2,scores,gbestIdx,vmax,wIter},vmax=vClampFrac*(ub-lb);particles=Transpose@Table[RandomReal[{lb[[i]],ub[[i]]},nParticles],{i,dim}];velocities=ConstantArray[0.,{nParticles,dim}];pbest=particles;pbestScore=lossFunc/@particles;gbestIdx=First@Ordering[pbestScore,1];gbestScore=pbestScore[[gbestIdx]];gbest=pbest[[gbestIdx]];Do[wIter=wMax-(wMax-wMin)*(j-1)/(nIter-1);r1=RandomReal[1,{nParticles,dim}];r2=RandomReal[1,{nParticles,dim}];velocities=wIter*velocities+c1*r1*(pbest-particles)+c2*r2*(ConstantArray[gbest,nParticles]-particles);velocities=clampMatrix[velocities,vmax];particles+=velocities;particles=clipMatrix[particles,lb,ub];scores=lossFunc/@particles;pbest=MapThread[If[#3<#2,#1,#4]&,{particles,pbestScore,scores,pbest}];pbestScore=MapThread[Min,{pbestScore,scores}];gbestIdx=First@Ordering[pbestScore,1];If[pbestScore[[gbestIdx]]<gbestScore,gbestScore=pbestScore[[gbestIdx]];gbest=pbest[[gbestIdx]];];AppendTo[convergence,gbestScore];If[j<10||Mod[j,10]==0,Print["PSO Iter ",j," Loss = ",gbestScore]];,{j,nIter}];{gbest,gbestScore,convergence}];
Dimension of search space
Bounds: P2,P3,P4 all strictly between P1 and Pn
Real gas case
Run PSO on WtotalSRK
Stage-by-stage detail
Convergence plot (log scale)
Ideal gas case
Run PSO on WtotalSRK
Stage-by-stage detail
Convergence plot (log scale)
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Optimal inter-stage pressures for a 4-stage propane refrigeration compressor train via PSO
by Housam Binous
Wolfram Community, STAFF PICKS, July 5, 2026
https://community.wolfram.com/groups/-/m/t/3744714
by Housam Binous
Wolfram Community, STAFF PICKS, July 5, 2026
https://community.wolfram.com/groups/-/m/t/3744714