Hello all,
I’m trying to rewrite by myself a very simple refinery model to improve my knowledge of GAMS. I’m using Marco oil refinery model found on the GAMS website as a starting point (https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_marco.html). I’m trying to make it a bit simpler (labels to help me understand what every equation does and a simpler process to start with). However, at one point I’m having a problem with an endogenous variable apparently and I don’t understand why. Here is my code (some equations are commented out):
option Limrow = 100;
Set component /SR-NAP, SR-GO, SR-FO/
crude /CrudeA, CrudeB /
specs /Density, Sulphur, VBN /
final /GASOLINE, GASOIL, FUELOIL /;
Positive Variable RunsMT(crude);
Parameter
Production(crude, component) '%WT of cut' / CrudeA.SR-NAP 0.2107
CrudeA.SR-GO 0.4820
CrudeA.SR-FO 0.3073
CrudeB.SR-NAP 0.2211
CrudeB.SR-GO 0.4217
CrudeB.SR-FO 0.3572 /
Price(final) 'Final product price in $/mt' /GASOLINE 580, GASOIL 550, FUELOIL 400 /
PriceCrude(crude) 'Crude price in $/bbl' /CrudeA 60, CrudeB 62/
CrudeAvails(crude) 'Max kbbls per crude' /CrudeA 300, CrudeB 300/;
Table Properties(crude, component, specs) 'Qualities of cuts'
Density Sulphur VBN
CrudeA.SR-NAP 0.712 0.0001
CrudeA.SR-GO 0.836 0.0001
CrudeA.SR-FO 0.860 0.5000 36
CrudeB.SR-NAP 0.600 0.0001
CrudeB.SR-GO 0.836 0.0001
CrudeB.SR-FO 0.858 0.5000 36;
Table SpecCrude(crude, specs) 'Qualities of crudes'
Density Sulphur
CrudeA 0.834414 0.1497
CrudeB 0.823392 0.0571;
Set limit 'Min or Max limit' /min, max/;
Table SpecLimits(limit, final, specs) 'Min Max Qualities for sales'
Density Sulphur VBN
max.GASOLINE 0.755 0.1
max.GASOIL 0.845 0.1
max.FUELOIL 5.991 3.5000 37
min.FUELOIL 0.800 32;
Set BlendingStream(final, component) 'blending possibility '
/ GASOLINE. (SR-NAP)
GASOIL. (SR-GO)
FUELOIL. (SR-GO, SR-FO)
/;
Parameter
BFComponent(crude, component) 'Component Barrel Factor'
BFCrude(crude) 'Crude Barrel Factor';
BFComponent(crude, component) = 1000 * (158.9872972 * Properties(crude, component, 'Density')) **(-1);
BFCrude(crude) = 1000 * (158.9872972 * SpecCrude(crude, 'Density')) **(-1);
display BFComponent, BFCrude;
Variable GPW, ComponentProperty, ComponentOrigin, BlendBalance, Blending, BlendingProperty, BlendingBF;
Positive variable Blend, SaleMT, SaleBBL, TotRunMT, TotRunBBL, ComponentAvail;
TotRunBBL.up = 300;
Equations OBJ,
MaxCrudePurchase,
SalesMT,
* SalesBBL,
TotRunsMT,
TotRunsBBL,
MaterialBalance,
ComponentAvails,
ComponentProperties,
ComponentOrigins
BlendingProcess
* BlendingProperties
BlendingBalance
* MaxSpecLimits
* BlendingBFs;
;
OBJ.. GPW =e= sum(final, SaleMT(final) * Price(final)) - sum(crude, RunsMT(crude) * PriceCrude(crude) * BFCrude(crude));
MaterialBalance.. sum(final, SaleMT(final)) - sum(crude, RunsMT(crude)) =e= 0;
BlendingBalance(component).. ComponentAvail(component) - sum(final, Blend(final, component)) =e= 0;
TotRunsMT.. TotRunMT =e= sum(crude, RunsMT(crude));
TotRunsBBL.. TotRunBBL =e= sum(crude, RunsMT(crude) * BFCrude(crude));
MaxCrudePurchase(crude).. RunsMT(crude) * BFCrude(crude) =l= CrudeAvails(crude);
ComponentAvails(component).. ComponentAvail(component) =e= sum(crude, ComponentOrigin(crude, component));
ComponentOrigins(crude, component).. ComponentOrigin(crude, component) =e= Production(crude, component) * RunsMT(crude);
ComponentProperties(component, specs).. ComponentProperty(component, specs) =e= ComponentAvail(component) * sum(crude, Properties(crude, component, specs));
BlendingProcess(final).. Blending(final) =e= sum(component$BlendingStream(final, component), Blend(final, component));
*BlendingProperties(final, specs).. BlendingProperty(final, specs) =e= sum(component, ComponentProperty(component, specs) * Blend(final, component));
*MaxSpecLimits(final, specs)$SpecLimits("max", final, specs).. BlendingProperty(final, specs) =l= SpecLimits("max", final, specs) * SaleMT(final);
SalesMT(final).. SaleMT(final) =e= Blending(final);
*BlendingBFs(final).. BlendingBF(final) =e= 1000 * (158.9872972 * BlendingProperty(final, 'Density')) ** (-1);
*SalesBBL(final).. SaleBBL(final) =e= SaleMT(final) * BlendingBF(final);
Model RefySimple /ALL/;
SOLVE RefySimple using LP maximizing GPW;
Parameter final_spec_component
final_spec_sales;
final_spec_component(component, specs) = ComponentProperty.l(component, specs) / ComponentAvail.l(component);
*final_spec_sales(final, specs) = BlendingProperty.l(final, specs) / SaleMT.l(final);
*display final_spec_component, final_spec_sales
The equation that is giving me the problem is the BlendingProperties one. Could someone tell me why is that? I am sure my model is linear: I take the property for every component and multiply it for the amount of that component used in the blend through the Blend variable. If I’m doing something wrong, could you point me in the right direction? The Marco model is doing everything with the LP model so I’m sure it’s possible to express this in a linear way. Thanks!