Hi,
I am trying to speed up scenario analysis by using GUSS (GAMS version 35.2.0), however I’m finding that many of the scenarios finish with model status 4 Infeasible, even though they are feasible (i.e., model status 1 Optimal) when solved in isolation. I’m banging my head against a brick wall with this and starting to think that I should have just set some loops running! It’s an LP problem, with CPLEX 20.1.0.1 as the solver.
Any advice would be much appreciated.
The code is pasted below and also attached, though it’s quite long. As I have no idea what might be the problem, I’ve been struggling to think of how to put together a shorter working example.
However, I have reduced the data in the csv files to just two scenarios and I’ve found that the order of the scenarios affects whether one or both of the scenarios solves, which seems quite strange to me.
The attached csv files are in the order where the first scenario solves and the second scenario is infeasible. Swapping the entries in deltaAvgAll.csv, minFurnaceFracAll.csv, and bVariableScrapChargeAll.csv has the effect of causing both scenarios to solve. (The two rows in windCFsAll.csv are identical as far as I know, similarly in solarCFsAll.csv.)
$ontext
main5.gms
Changes from v4:
- In correctly-named project.
Changes from v3:
- Loops replaced with GUSS approach to scenario analysis, in an attempt to speed up solution.
Optimisation tool for hydrogen-based steelmaking (H-DRI+EAF),
including wind and solar power, electricity storage, electrolysers, and hydrogen storage.
Assumptions/limitations:
- Time intervals are all of equal length (i.e., tInt is a scalar)
- tInt is currently hardcoded as 1 hour (rather than using timestamps in renewables data)
- AC/DC conversion is ignored in battery storage / solar calculations
- There is no limit to maximum battery C rate
- No material storage is included
- There is no lower limit on shaft furnace and EAF turn-down
- EAF can turn-down continuously (rather than just running / not running)
- No heat loss incorporated when turning down
- Perfect foresight of renewables availabilities
Development ideas:
- Shaft furnace cost could be bundled in with EAF cost, assuming both can vary together.
- Could do a parameter sweep or wind/solar fraction to see what effect they make on costs and savings from flexibility.
$offtext
Set components 'Components of the system'
/ EAF 'in tCS/h'
shaftFurnace 'in tDRI/h'
electrolyser 'in MW.H2'
fuelCell 'in MW'
wind 'in MW'
solar 'in MW'
batteryStorageEnergy 'in MWh'
batteryStoragePower 'in MW'
hydrogenStorageEnergy 'in MWh'
hydrogenStoragePower 'in MW' /
t 'Time steps for wind and solar data';
* Efficiencies need checking.
* For now, could just set batteryStoragePower component in netcosts to 0 GBP and limit charge/discharge power by maximum C-rates (e.g., 0.3 and 1).
* Time interval tInt is very important and so should be set in a different way.
* EAF and electrolyser costs from: V. Vogl, M. Åhman, L.J. Nilsson, Assessment of hydrogen direct reduction for fossil-free steelmaking, J. Clean. Prod., 203 (2018), pp. 736-745
* Shaft furnace and fuel cell costs from: A. Toktarova et al, Interaction between electrified steel production and the north European electricity system, https://www.sciencedirect.com/science/article/pii/S0306261922000654
* Appendix of this paper suggests shaft furnace costs are per tonne of DRI (not crude steel): M. Xylia et al, Weighing regional scrap availability in global pathways for steel production processes, https://link.springer.com/article/10.1007/s12053-017-9583-7?utm_source=getftr&utm_medium=getftr&utm_campaign=getftr_pilot
* Wind and solar costs from: BEIS Electricity Generation Costs 2020
* Battery costs per kWh from: Lazard's Levelized Cost of Storage (https://www.lazard.com/media/451882/lazards-levelized-cost-of-storage-version-70-vf.pdf) and 2025 Moderate estimate in NREL 2022 Annual Technology Baseline
* Battery costs per kW from: 2025 Moderate estimate from NREL 2022 Annual Technology Baseline
Parameter netcosts(components) 'Costs in GBP per unit installed capacity - THESE NEED CHECKING'
/ EAF 163, shaftFurnace 284, electrolyser 522e3, fuelCell 397, wind 3104e3, solar 617e3
batteryStorageEnergy 179e3, batteryStoragePower 198e3, hydrogenStorageEnergy 670, hydrogenStoragePower 0 /
windCapacityFactors(t) 'Wind capacity factors, 0 to 1'
solarCapacityFactors(t) 'Solar capacity factors, 0 to 1'
avgScrapChargeFraction 'Scrap charge fractions' / 0 /
minFurnaceTurnDownFrac 'Furnace turn-down scenarios' / 1 /
bVariableScrapCharge 'Scrap charge variability scenarios' / 0 /;
* avgScrapChargeFraction 'Scrap charge fractions' / 0.5 /
* minFurnaceTurnDownFrac 'Furnace turn-down scenarios' / 1 /
* bVariableScrapCharge 'Scrap charge variability scenarios' / 0 /;
* Convert furnace costs from £/tCS/yr to £/tCS/h.
netcosts('EAF') = netcosts('EAF') * 8766;
netcosts('shaftFurnace') = netcosts('shaftFurnace') * 8766;
Scalars electrolyserEfficiency /0.72/
fuelCellEfficiency /0.6/
batteryChargeEfficiency /0.92/
batteryDischargeEfficiency /0.92/
hydrogenStorageChargeEfficiency /1/
hydrogenStorageDischargeEfficiency /1/
annualOutput_MtCS 'in MtCS/yr' / 1 /
tInt 'in hours' / 3 /
r 'Discount rate' / 0.08 /
n 'Project life in years' / 20 /
CRF 'Capital recovery factor';
* Calculate capital recovery factor.
CRF = ( r * (1+r) ** n ) / ( (1+r) ** n - 1 );
* Load wind capacity factors.
*$call csv2gdx windCapacityFactors.csv id=windCapacityFactors index=1 values=2..lastCol useHeader=y trace=0
*$ifE errorLevel<>0 $abort Problems reading wind capacity factors file!
*$gdxIn windCapacityFactors.gdx
*$load t = dim1
*$load windCapacityFactors
*$gdxIn
*
** Load solar capacity factors.
*$call csv2gdx solarCapacityFactors.csv id=solarCapacityFactors index=1 values=2..lastCol useHeader=y trace=0
*$ifE errorLevel<>0 $abort Problems reading solar capacity factors file!
*$gdxIn solarCapacityFactors.gdx
** $load t = dim1
*$load solarCapacityFactors
*$gdxIn
* Load capacity factors and other parameters for the scenarios.
* Table allData(t,*) 'Wind/solar capacity factors, 0 to 1; Electricity demand of the site in MW';
Set scenariostorun 'scenarios';
Parameter windCFsAll(scenariostorun,t) 'Wind capacity factors for each scenario, 0 to 1';
Parameter solarCFsAll(scenariostorun,t) 'Solar capacity factors for each scenario, 0 to 1';
Parameter deltaAvgAll(scenariostorun) 'Average share of scrap in EAF charge for each scenario, 0 to 1';
Parameter minFurnaceFracAll(scenariostorun) 'Minimum furnace output for each scenario, 0 to 1';
Parameter bVariableScrapChargeAll(scenariostorun) 'Variable scrap charge flag for each scenario, 0 or 1';
$call csv2gdx windCFsAll.csv id=windCFsAll index=1 values=2..lastCol useHeader=y trace=0
$ifE errorLevel<>0 $abort Problems reading data file!
$gdxIn windCFsAll.gdx
$load scenariostorun = dim1
$load t = dim2
$load windCFsAll
$gdxIn
$call csv2gdx solarCFsAll.csv id=solarCFsAll index=1 values=2..lastCol useHeader=y trace=0
$ifE errorLevel<>0 $abort Problems reading data file!
$gdxIn solarCFsAll.gdx
$load solarCFsAll
$gdxIn
$call csv2gdx deltaAvgAll.csv id=deltaAvgAll index=1 value=2 useHeader=y trace=0
$ifE errorLevel<>0 $abort Problems reading data file!
$gdxIn deltaAvgAll.gdx
$load deltaAvgAll
$gdxIn
$call csv2gdx minFurnaceFracAll.csv id=minFurnaceFracAll index=1 value=2 useHeader=y trace=0
$ifE errorLevel<>0 $abort Problems reading data file!
$gdxIn minFurnaceFracAll.gdx
$load minFurnaceFracAll
$gdxIn
$call csv2gdx bVariableScrapChargeAll.csv id=bVariableScrapChargeAll index=1 value=2 useHeader=y trace=0
$ifE errorLevel<>0 $abort Problems reading data file!
$gdxIn bVariableScrapChargeAll.gdx
$load bVariableScrapChargeAll
$gdxIn
Display scenariostorun, t, deltaAvgAll, minFurnaceFracAll, bVariableScrapChargeAll, windCFsAll, solarCFsAll;
* Set capacity factors for base case.
windCapacityFactors(t) = windCFsAll('s0',t);
solarCapacityFactors(t) = solarCFsAll('s0',t);
Positive Variables capacities(components) 'Unit capacities'
steelOutput(t) 'in tCS/h'
DRI(t) 'in tDRI/h'
scrap(t) 'in t/h'
feInScrap(t) 'in t/h'
feInDRI(t) 'in t/h'
H2fromElectrolyser(t) 'in MW.H2'
powerToEAF(t) 'in MW'
powerToHydrogenHeater(t) 'in MW'
windOutput(t) 'in MW'
solarOutput(t) 'in MW'
powerToBattery(t) 'in MW'
powerFromBattery(t) 'in MW'
powerToElectrolyser(t) 'in MW'
powerFromFuelCell(t) 'in MW'
powerToHydrogenStorage(t) 'in MW'
powerFromHydrogenStorage(t) 'in MW'
E0battery 'in MWh'
E0hydrogenStorage 'in MWh'
energyInBattery(t) 'in MWh'
energyInHydrogenStorage(t) 'in MWh';
Variables cost 'Total sum of costs'
annualCost 'Annualised total costs'
costPerTCS 'Cost per tonne of crude steel';
Equations costAcct 'Cost accounting equation (to be minimised)'
annualCostAcct 'Annual cost accounting equation'
costPerTCSAcct 'Cost per tCS accounting equation'
powerToEAFEqn(t) 'Power to EAF'
powerToHydrogenHeaterEqn(t) 'Power to hydrogen heater'
electricityBalanceEqn(t) 'Balance electricity flows'
hydrogenBalanceEqn(t) 'Balance H2 flows'
batteryChargeLimitEqn(t) 'Limit storage powers'
batteryDischargeLimitEqn(t)
hydrogenStorageChargeLimitEqn(t)
hydrogenStorageDischargeLimitEqn(t)
windEqn(t) 'Limit wind output by wind capacity factors'
solarEqn(t) 'Limit solar output by solar capacity factors'
batteryEnergyEqn(t) 'Energy in battery at each timestep'
hydrogenStorageEnergyEqn(t) 'Energy in hydrogen storage at each timestep'
batteryEnergyUpperLimit(t) 'Limit energy in battery storage'
hydrogenStorageEnergyUpperLimit(t) 'Limit energy in hydrogen storage'
E0batteryLimit 'Initial energy in battery less than capacity'
E0hydrogenStorageLimit 'Initial energy in hydrogen storage less than capacity'
terminalConstraintBattery 'Constraint on energy in battery at end of period'
terminalConstraintHydrogenStorage 'Constraint on energy in hydrogen storage at end of period'
steelProductionEqn 'Ensure steel production requirement is met'
feInScrapEqn(t) 'Calculate iron content in scrap'
feInDRIEqn(t) 'Calculate iron content in DRI'
ironBalanceEqn(t) 'Ensure iron flow is sufficient for the flow of steel'
EAFcapacityUBEqn(t) 'Limit steel output by EAF capacity'
shaftFurnaceUBEqn(t) 'Limit DRI output by shaft furnace capacity'
EAFcapacityLBEqn(t)
shaftFurnaceLBEqn(t)
scrapChargeFractionUBEqn(t) 'Limit scrap charge to be less than the required average (if no flexibility) or 1'
scrapChargeAverageEqn 'Ensure average scrap charge meets the requirement'
electrolyserEqn(t) 'Limit electrolyser power by electrolyser capacity'
fuelCellEqn(t) 'Limit fuel cell output by fuel cell capacity';
* Define objective function (total cost) and additional cost calculations.
costAcct.. cost =e= sum(components, netcosts(components)*capacities(components));
annualCostAcct.. annualCost =e= cost * CRF;
costPerTCSAcct.. costPerTCS =e= annualCost / ( annualOutput_MtCS * 1e6 );
* Limit flows to/from storage.
batteryChargeLimitEqn(t).. powerToBattery(t) =l= capacities('batteryStoragePower');
batteryDischargeLimitEqn(t).. powerFromBattery(t) =l= capacities('batteryStoragePower');
hydrogenStorageChargeLimitEqn(t).. powerToHydrogenStorage(t) =l= capacities('hydrogenStoragePower');
hydrogenStorageDischargeLimitEqn(t).. powerFromHydrogenStorage(t) =l= capacities('hydrogenStoragePower');
* Limit energy in storage.
alias (t, tAlias);
* energyInBattery(t) .. E0battery + sum(tAlias$[ord(tAlias) * tInt; * Could potentially simplify in this way.
batteryEnergyEqn(t).. energyInBattery(t) =e= E0battery + sum(tAlias$(ord(tAlias) le ord(t)), tInt * (powerToBattery(tAlias) * batteryChargeEfficiency - powerFromBattery(tAlias) / batteryDischargeEfficiency));
hydrogenStorageEnergyEqn(t).. energyInHydrogenStorage(t) =e= E0hydrogenStorage + sum(tAlias$(ord(tAlias) le ord(t)), tInt * (powerToHydrogenStorage(tAlias) * hydrogenStorageChargeEfficiency - powerFromHydrogenStorage(tAlias) / hydrogenStorageDischargeEfficiency));
batteryEnergyUpperLimit(t).. energyInBattery(t) =l= capacities('batteryStorageEnergy');
hydrogenStorageEnergyUpperLimit(t).. energyInHydrogenStorage(t) =l= capacities('hydrogenStorageEnergy');
* Limit initial energy in storage.
E0batteryLimit.. E0battery =l= capacities('batteryStorageEnergy');
E0hydrogenStorageLimit.. E0hydrogenStorage =l= capacities('hydrogenStorageEnergy');
* Terminal constraints on energy in storage, ensuring that SOC at end of analysis period == SOC at start.
terminalConstraintBattery.. sum(t, energyInBattery(t) - E0battery) =e= 0;
terminalConstraintHydrogenStorage.. sum(t, energyInHydrogenStorage(t) - E0hydrogenStorage) =e= 0;
* Ensure electricity demand is met.
* Values generally from Vogl (2018) and Devlin papers (2022 & 2023)
* Power to EAF should be match up with Vogl numbers at the extremes of scrap charge ratio, i.e., 2.4 GJ/tLS (667 kWh/tLS) when delta = 1 and 2.711 GJ/tLS (753 kWh/tLS) when delta = 0
* 1606.1 MJ/tDRI hydrogen heating requirement from Devlin & Yang (2022)
*powerToEAFEqn(t).. powerToEAF(t) =e= steelOutput(t) * ( 2.4 + 0.85 * 3.6 * (1 - scrapChargeFraction(t)) / 10 ) / 3.6 / tInt;
powerToEAFEqn(t).. powerToEAF(t) =g= ( 2.4 * feInScrap(t) + 2.711 * feInDRI(t) ) / 3.6;
powerToHydrogenHeaterEqn(t).. powerToHydrogenHeater(t) =g= 1606.1 * DRI(t) / 3.6e3;
electricityBalanceEqn(t)..
windOutput(t) + solarOutput(t) + powerFromBattery(t) - powerToBattery(t) - powerToElectrolyser(t) + powerFromFuelCell(t) =g= powerToEAF(t) + powerToHydrogenHeater(t);
*electricityBalanceEqn(t)..
* windOutput(t) + solarOutput(t) + powerFromBattery(t) - powerToBattery(t) - powerToElectrolyser(t) + powerFromFuelCell(t) =g= 1.1 * steelOutput(t) / tInt;
* Ensure hydrogen demand is met.
hydrogenBalanceEqn(t)..
powerToElectrolyser(t) * electrolyserEfficiency - powerFromFuelCell(t) / fuelCellEfficiency + powerFromHydrogenStorage(t) - powerToHydrogenStorage(t) =g= 2.4 * electrolyserEfficiency * feInDRI(t);
* Ensure annual steel production requirement is met, converting RHS from MtCS to tCS.
steelProductionEqn..
sum(t, steelOutput(t)*tInt) =g= annualOutput_MtCS * 1e6 * sum(t, tInt) / 8766;
* Limit scrap charge fraction to either avgScrapChargeFraction (if bVariableScrapCharge is false) or 1 (if bVariableScrapCharge is true).
scrapChargeFractionUBEqn(t).. scrap(t) =l= (scrap(t) + DRI(t)) * (avgScrapChargeFraction + (1-avgScrapChargeFraction) * bVariableScrapCharge);
scrapChargeAverageEqn.. sum(t, scrap(t)*tInt) * (1 - avgScrapChargeFraction) =l= sum(t, DRI(t)*tInt) * avgScrapChargeFraction;
* Ensure iron flow is sufficient for the flow of steel (assuming no DRI storage and that steel is 100% iron).
* Coefficients calculated using voglModelCalcs.ipynb in Google Colab (file in Google Drive)
feInScrapEqn(t).. feInScrap(t) =l= 0.95 * scrap(t);
feInDRIEqn(t).. feInDRI(t) =l= 0.915 * DRI(t);
ironBalanceEqn(t).. feInScrap(t) + feInDRI(t) =g= steelOutput(t);
* Ensure all energy and resource flows are within unit capacities.
EAFcapacityUBEqn(t)..
steelOutput(t) =l= capacities('EAF');
shaftFurnaceUBEqn(t)..
DRI(t) =l= capacities('shaftFurnace');
EAFcapacityLBEqn(t)..
steelOutput(t) =g= minFurnaceTurnDownFrac * capacities('EAF');
shaftFurnaceLBEqn(t)..
DRI(t) =g= minFurnaceTurnDownFrac * capacities('shaftFurnace');
windEqn(t).. windOutput(t) =l= windCapacityFactors(t) * capacities('wind');
solarEqn(t).. solarOutput(t) =l= solarCapacityFactors(t) * capacities('solar');
electrolyserEqn(t).. powerToElectrolyser(t) =l= capacities('electrolyser');
fuelCellEqn(t).. powerFromFuelCell(t) =l= capacities('fuelCell');
* Set up the models and select which constraint equations to include.
Model steelProblem /all/;
$ontext
* Prepare a csv file to save the results.
File results / resultsFile.csv /;
results.pc = 5;
put results;
put 'minFurnaceFrac', 'scrapChargeFrac', 'variableScrapCharge', 'cost', 'costPerTCS';
loop(components, put components.tl);
put /;
* Step through each location/year.
loop(locationsYears,
currfile = inputfiles(locationsYears);
*$call csv2gdx dataJan2015.csv id=allData index=1 values=2..lastCol useHeader=y trace=0
$call csv2gdx currfile id=allData index=1 values=2..lastCol useHeader=y trace=0
$ifE errorLevel<>0 $abort Problems reading data file!
*$gdxIn dataJan2015.gdx
$gdxIn currfile
$load t = dim1
$load allData
$gdxIn
* Solve the model and save the results.
solve steelProblem using LP minimizing cost;
put minFurnaceTurnDownFrac, avgScrapChargeFraction, bVariableScrapCharge, cost.l, costPerTCS.l;
loop(components, put capacities.l(components));
put /;
);
putclose;
$offtext
Parameter resultantCostPerTCS(scenariostorun) 'Collector for cost per tonne of crude steel'
scopt / SkipBaseCase 1 /;
* Define the scenarios to be analysed.
set dict / scenariostorun .scenario .''
windCapacityFactors .param .windCFsAll
solarCapacityFactors .param .solarCFsAll
avgScrapChargeFraction .param .deltaAvgAll
minFurnaceTurnDownFrac .param .minFurnaceFracAll
bVariableScrapCharge .param .bVariableScrapChargeAll
costPerTCS .level .resultantCostPerTCS
scopt .opt .''
/;
* gussoptions .opt .solutionstatus
* Solve the model for all scenarios.
solve steelProblem using LP minimizing cost scenario dict;
Thanks,
Andrew
solarCFsAll.csv (76.8 KB)
windCFsAll.csv (101 KB)
main5.gms (18.7 KB)
bVariableScrapChargeAll.csv (29 Bytes)
deltaAvgAll.csv (24 Bytes)
minFurnaceFracAll.csv (32 Bytes)