Closed dguittet closed 1 year ago
Hi @adowling2, I tried the Degeneracy Hunter but at this step: dh = DegeneracyHunter(m, solver=pyo.SolverFactory('cbc'))
I run into this error:
if not AmplInterface.available():
raise RuntimeError(
"Cannot load the PyNumero ASL interface (pynumero_ASL)")
@dguittet Are you running this on macOS, Windows, or Linux? If Windows or Linux, you want to install the IDAES extensions which include the "good" Ipopt (with HSL linear algebra). This will also give you the pynumero_ASL library. For macOS users, it is much more complicated. I recommend VirtualBox + Ubuntu.
Also, using "good" Ipopt may dramatically speed up the solution.
@xiangao1 Can you share with @dguittet the code we have for initializing an IDAES model from a library of JSONs?
@adowling2 Great, I was able to install Pynumero from Pyomo and get dh
to work! I'm on Mac but I built IPOPT with HSL. However, the timing results I showed last week are from Linux using idaes-ext (HPC not locally).
So after running IPOPT with max_its = 0, there does show up constraint violations. This is the code here:
# no infeasible constraints or bounds show up during this block
solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties")
log_infeasible_constraints(m, logger=solve_log, tol=1e-4, log_expression=True, log_variables=True)
log_infeasible_bounds(m, logger=solve_log, tol=1e-4)
time_to_create_model = default_timer() - start
opt.options['max_iter'] = 0
status_obj, solved, iters, time, regu = ipopt_solve_with_stats(m, opt, opt.options['max_iter'], 60*210)
# Now violations do show up, replicates what dh.check_residuals(tol=0.1) shows below
log_infeasible_constraints(m, logger=solve_log, tol=1e-4, log_expression=True, log_variables=True)
log_infeasible_bounds(m, logger=solve_log, tol=1e-4)
dh = DegeneracyHunter(m, solver=pyo.SolverFactory('cbc'))
dh.check_residuals(tol=0.1)
dh.check_variable_bounds(tol=1.0)
Is there anything that could or should be done here? I tried re-initializing the model, but that gives a pyomo error:
After the full solve, there is some small residuals, within the tolerance set for IPOPT and check_rank_equality_constraints
gives:
Checking rank of Jacobian of equality constraints...
Model contains 537 equality constraints and 556 variables.
Computing the 10 smallest singular value(s)
Since @adowling2 is away, @jsiirola do you have any suggestions?
Here's a slightly more detailed time profiling of solving for 4380 hrs at a time. Total run time was 6.25 hours: IPOPT solve was 5.27 hours, pre-solve was 0.829 hours. Within pre-solve, initialization was 0.55 hours and cloning was 0.269 hours. IPOPT time (in the opt.solve
call) was 84% of the total time.
{
"h2_price": 2.0,
"ts": 4380,
"wind_mw": 200.0,
"batt_mw": 330.598,
"pem_mw": 0.000441,
"tank_kgH2": 0.0737,
"turb_mw": 0.000208,
"avg_turb_eff": 0.329,
"annual_rev_h2": -8.467,
"annual_rev_E": 50944289,
"NPV": 388573302,
"idaes_creation_time": 970,
"idaes_init_time": 2010,
"multiperiod_init": 0.000006,
"multiperiod_build": 2982,
"model_vars": 1213262,
"ipopt_solved": "True",
"ipopt_its": 2368,
"ipopt_time": 19001,
"ipopt_reg": "-19.1",
"presolve_time": 2986,
"solve_time": 19001,
"total_run_time": 22525
}
I've uploaded the log from IPOPT as well. Overall, there aren't restoration phases and there's some regularization at the end. 4380_hr.log.zip The time split is:
Total CPU secs in IPOPT (w/o function evaluations) = 15316.738
Total CPU secs in NLP function evaluations = 3685.171
At $2.0/kg H2, there's no PEM selected, so this has finished faster than the other runs which have higher H2 prices and more PEM for the optimal design. I also have 1 year runs going right now. I'll update the rest as they finish.
And here's the time profiling for 1 one of the 8736 hr runs (the other 2 are still incomplete). Total run time was 21.76 hrs, 19.02 hrs of which was in IPOPT, or 87%.
{
"ts": 8736,
"h2_price": 3.0,
"wind_mw": 200.0,
"batt_mw": 118.4388,
"pem_mw": 37.9098,
"tank_kgH2": 0.1515,
"turb_mw": 0.0001835,
"avg_turb_eff": 0.3253,
"annual_rev_h2": 12095933,
"annual_rev_E": 17860958,
"NPV": 257857205,
"idaes_creation_time": 1803,
"idaes_init_time": 6816,
"multiperiod_init": 0.0000056,
"multiperiod_build": 8622,
"model_vars": 2419874,
"ipopt_solved": "True",
"ipopt_its": 2607,
"ipopt_time": 68485.542,
"ipopt_reg": "-",
"presolve_time": 8629,
"solve_time": 68485,
"total_run_time": 78347
}
@dguittet any news on this?
@dguittet I just had a chance to look at the logs:
2092 -2.5423391e+03 1.58e-01 6.31e+00 -7.0 5.13e+01 -7.2 4.31e-02 8.76e-01h 1
2093 -2.5423662e+03 2.43e-01 4.61e+00 -7.0 5.64e+02 - 3.03e-01 6.49e-01h 1
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2094 -2.5423870e+03 1.50e-01 4.65e+00 -7.0 1.79e+02 - 6.43e-02 3.77e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2095 -2.5424095e+03 7.93e-02 4.46e+00 -7.0 1.90e+02 - 8.82e-02 3.73e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2096 -2.5424355e+03 4.55e-02 3.64e+00 -7.0 2.67e+02 - 2.06e-01 3.78e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2097 -2.5424683e+03 5.32e-02 3.64e+00 -7.0 2.97e+02 - 3.48e-02 3.62e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2098 -2.5425033e+03 5.89e-02 2.77e+00 -7.0 2.81e+02 - 2.44e-01 3.60e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2099 -2.5425477e+03 5.99e-02 2.42e+00 -7.0 2.60e+02 - 1.28e-01 3.35e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
2100 -2.5425977e+03 5.59e-02 1.97e+00 -7.0 2.40e+02 - 1.76e-01 3.23e-01h 2
2101 -2.5427031e+03 5.98e-02 1.98e+00 -7.0 3.05e+02 - 6.43e-03 5.60e-01h 1
2102 -2.5427342e+03 5.09e-02 1.81e+00 -7.0 3.32e+02 - 9.70e-02 1.63e-01h 1
2103 -2.5427676e+03 4.25e-02 1.77e+00 -7.0 3.76e+02 - 4.78e-02 1.59e-01f 1
2104 -2.5427993e+03 3.83e-02 1.48e+00 -7.0 4.03e+02 - 1.64e-01 1.44e-01f 1
2105 -2.5428366e+03 3.63e-02 1.10e+00 -7.0 4.90e+02 - 2.32e-01 1.44e-01h 1
2106 -2.5428738e+03 3.45e+00 1.11e+00 -7.0 1.04e+04 - 2.34e-02 1.15e-01h 1
2107 -2.5429087e+03 3.10e+00 9.24e-01 -7.0 1.38e+03 - 1.51e-01 1.03e-01h 1
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2108 -2.5430720e+03 1.65e+00 8.80e-01 -7.0 4.01e+01 -7.6 1.52e-01 4.69e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2109 -2.5431932e+03 1.24e+00 7.01e-01 -7.0 1.01e+03 - 2.16e-01 2.66e-01h 2
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
My suggestion is to enable the halt_on_ampl_error
option. Then, set symbolic_solver_labels=True
in the solve
call and rerun. This will tell us which constraint is causing the ampl error. This normally points to a model problem, that when sometimes addressed with a reformulation, leads to much faster Ipopt convergence.
@adowling2 Yes, I've tracked that problem down using those options and it comes down to the isentropic_properties block in the h2 turbine's compressor and turbine, where there are log functions of mole fraction ratios. Since H2 is not optimal, as it gets minimized, IPOPT will try out values where there is a very, small negative number for those mole fraction ratios of hydrogen, leading to the ampl_error.
At @bknueven's suggestion, I have tried running IPOPT with bound_relax_factor= 0.0
, and this prevents those AMPL errors but results in practically the same number of iterations and solve time.
@dguittet that makes sense. @radhakrishnatg does it make sense to bound these mole fractions inside the H2 turbine model components? IMO, this sounds like something better to fix in the model than on the solver side.
@ksbeattie This could be a low-priority issue. It is not needed for a milestone this year.
@adowling2 @radhakrishnatg I've tried bounding those mole fractions as well but it didn't improve the situation. I think it may have taken longer or had other solve problems.
EDIT: Though I haven't tried bounding the mole fractions since fixing the air feed ratio which had been a big issue before. It's possible that now that the air feed ratio is fixed, other bounds could be tightened and have a positive impact.
@adowling2 The plan is to make a small table of optimal designs by front loading the simulations. Then the double loop code can just read the optimal sizes given a particular H2 price (or some other parameter).
@ksbeattie This could be a low-priority issue. It is not needed for a milestone this year.
@adowling2, thanks. I've moved this to the June 2022 release board, for now.
@dguittet that makes sense. @radhakrishnatg does it make sense to bound these mole fractions inside the H2 turbine model components? IMO, this sounds like something better to fix in the model than on the solver side.
@adowling2 @dguittet This is how I initialize the translator, mixer and hydrogen turbine models:
m.fs.translator.outlet.mole_frac_comp[0, "hydrogen"].fix(0.99)
m.fs.translator.outlet.mole_frac_comp[0, "oxygen"].fix(0.01 / 4)
m.fs.translator.outlet.mole_frac_comp[0, "argon"].fix(0.01 / 4)
m.fs.translator.outlet.mole_frac_comp[0, "nitrogen"].fix(0.01 / 4)
m.fs.translator.outlet.mole_frac_comp[0, "water"].fix(0.01 / 4)
propagate_state(m.fs.arc_h2_tank_to_translator)
m.fs.translator.initialize()
# Fix the degrees of freedom of mixer and initialize
m.fs.mixer.air_feed.flow_mol[0].fix(
m.fs.h2_tank.outlet_to_turbine.flow_mol[0].value * air_h2_ratio
)
m.fs.mixer.air_feed.temperature[0].fix(pem_outlet_temperature)
m.fs.mixer.air_feed.pressure[0].fix(pem_outlet_pressure * 1e5)
m.fs.mixer.air_feed.mole_frac_comp[0, "oxygen"].fix(0.2054)
m.fs.mixer.air_feed.mole_frac_comp[0, "argon"].fix(0.0032)
m.fs.mixer.air_feed.mole_frac_comp[0, "nitrogen"].fix(0.7672)
m.fs.mixer.air_feed.mole_frac_comp[0, "water"].fix(0.0240)
m.fs.mixer.air_feed.mole_frac_comp[0, "hydrogen"].fix(2e-4)
propagate_state(m.fs.arc_translator_to_mixer)
m.fs.mixer.initialize()
After initializing the model, I unfix the molar flowrate of air to the turbine m.fs.mixer.air_feed.flow_mol.unfix()
and add a constraint fixing the air-to-fuel ratio:
# Add a constraint to maintain the air to hydrogen flow ratio
m.fs.mixer.air_h2_ratio = Constraint(
expr=m.fs.mixer.air_feed.flow_mol[0] ==
air_h2_ratio * m.fs.mixer.hydrogen_feed.flow_mol[0])
However, I do not unfix the mole fractions, because the composition of both the hydrogen feed and the air feed does not change during the operation. Because of this, the mole fraction of all the components entering the hydrogen turbine remain fixed, so I did not add any bounds on them. Then the only thing we need to worry about is the mole fraction of hydrogen in the outlet.
@dguittet are you fixing the conversion of hydrogen in the stoichiometric reactor? If not, please fix it to 0.99 m.fs.h2_turbine.stoic_reactor.conversion.fix(0.99)
. Otherwise, the optimizer will make it 0.999999 (basically, push it as close to 1 as possible), and the mole fraction of hydrogen starts approaching 0. There are two approaches to fix this. Either bound the mole fraction of hydrogen as @adowling2 suggested, or fix/bound the conversion of hydrogen. In the optimal solution, the conversion of hydrogen will always be 1, so we can fix it to a value close to 1.
Please let me know if this does not help. I can try a few more things.
I have an approach to reduce the model build time. I am not sure if it works, but we can definitely try it and see when you have some time.
@radhakrishnatg I also have the air feed ratio fixed, and the stoic_reactor's conversion rate to 0.99. I can re-run the model to report the AMPL errors and see if those are different now compared to when I last checked them, which was before the air feed ratio was fixed. If it's not mole fractions, it may be just the flow rates, perhaps.
The model build time is ~15% of the total time, most of it is around re-initializing the wind and wind_to_grid blocks to adapt to the time-varying wind resource. The rest of the unit models are not re-initialized.
Well, with the latest code, I can not reproduce that AMPL error even over a few different H2 prices and time horizons...
Okay. Is IPOPT still taking a long time to solve the problem? With the mole fraction and conversion values fixed, the hydrogen turbine model becomes linear. So, I think it may not be the reason for the slow convergence of IPOPT. Do you observe a significant reduction in computational time when you remove the hydrogen turbine?
Need to re-visit this soon to be able to annual simulations in a reasonable time, however, this doesn't need to be in the October 2022 release
@dguittet Could you please point me to the final version of the code? I would like to try a couple of things and see if they work.
Couldn't speed up the nonlinear flowsheet for annual runs. Instead linearized model here: https://github.com/dguittet/dispatches/blob/re_h2_storage/dispatches/case_studies/renewables_h2_case/wind_battery_hydrogen_flowsheet.py
Closing this issue for now and future PR will bring in the above file
Currently, the full flowsheet (Wind + Battery + PEM + H2 Tank + Turbine) has not beed solved for a full year of data. Speed up model to make this possible.
Difficulties with this model include:
Things implemented or tried already:
purchased_hydrogen_feed_state
existslog_infeasible_constraints
&log_infeasible_bounds
to reduce initial iteration errorThings to try: