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/ ... 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):
Code: Select all
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