Model status Infeasible when using GUSS

questions about GAMS' tools
Post Reply
Posts: 1
Joined: 3 months ago

Model status Infeasible when using GUSS

Post by andrewpimm »


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 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.)

Code: Select all

  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.
  - 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.
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,
* 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,
* Wind and solar costs from: BEIS Electricity Generation Costs 2020
* Battery costs per kWh from: Lazard's Levelized Cost of Storage ( 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
** 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

* 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

$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

$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

$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

$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

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'
                     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'
                     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;
      windOutput(t) + solarOutput(t) + powerFromBattery(t) - powerToBattery(t) - powerToElectrolyser(t) + powerFromFuelCell(t) =g= powerToEAF(t) + powerToHydrogenHeater(t);
*      windOutput(t) + solarOutput(t) + powerFromBattery(t) - powerToBattery(t) - powerToElectrolyser(t) + powerFromFuelCell(t) =g= 1.1 * steelOutput(t) / tInt;

* Ensure hydrogen demand is met.
      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.
      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.
    steelOutput(t) =l= capacities('EAF');

    DRI(t) =l= capacities('shaftFurnace');

    steelOutput(t) =g= minFurnaceTurnDownFrac * capacities('EAF');
    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/;

* 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;
put /;

* Step through each location/year.
    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
* 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 /;


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;
(18.74 KiB) Downloaded 18 times
(101.27 KiB) Downloaded 20 times
(76.81 KiB) Downloaded 18 times
(32 Bytes) Downloaded 18 times
(24 Bytes) Downloaded 17 times
(29 Bytes) Downloaded 19 times
User avatar
Posts: 1017
Joined: 7 years ago

Re: Model status Infeasible when using GUSS

Post by bussieck »

Not clear what problem you describe. When I run the model you uploaded I get:

Code: Select all

--- Scenario s0:        Records read: 4650
--- Solution statistics for scenario s0
---   Unmatched Records  0
---   Solver Status      1 Normal Completion
---   Model Status       1 Optimal
---   Objective Value    2.80346e+09
---   Resource Usage     281.796
--- Scenario s1:        Records read: 4651
--- Solution statistics for scenario s1
---   Unmatched Records  0
---   Solver Status      1 Normal Completion
---   Model Status       1 Optimal
---   Objective Value    1.39903e+09
---   Resource Usage     683.375
with an explicit loop like this

Code: Select all

   option clear=windCapacityFactors, clear=solarCapacityFactors;
   windCapacityFactors(t)  = windCFsAll(scenariostorun,t);
   solarCapacityFactors(t) = solarCFsAll(scenariostorun,t);
   avgScrapChargeFraction  = deltaAvgAll(scenariostorun);
   minFurnaceTurnDownFrac  = minFurnaceFracAll(scenariostorun);
   bVariableScrapCharge    = bVariableScrapChargeAll(scenariostorun);
   solve steelProblem using LP minimizing cost;
I get

Code: Select all

LOOPS                           scenariostorun   s0

**** SOLVER STATUS     1 Normal Completion
**** MODEL STATUS      1 Optimal
**** OBJECTIVE VALUE       2803458327.7206

 RESOURCE USAGE, LIMIT        277.672 10000000000.000

LOOPS                           scenariostorun   s1

**** SOLVER STATUS     1 Normal Completion
**** MODEL STATUS      1 Optimal
**** OBJECTIVE VALUE       1399028162.3276

 RESOURCE USAGE, LIMIT        647.234 10000000000.000
So I see no difference in execution. Can you provide a smoking gun?

Post Reply