idaes-pse
idaes-pse copied to clipboard
IDAES
Demonstrate PyROS on rankine cycle IES example
Overall goal: solve 24 hour bidding problem for Rankine cycle with uncertainty in LMPs using Pyros
Main Steps:
- [ ] Debug separation problem
- [ ] Specify uncertainty sets for 24-hour example
Code: https://github.com/shermanjasonaf/idaes-pse/tree/multiperiod-rankine-uncertainty/idaes/apps/multiperiod/examples
-
multiperiod_rankine_uncertainty.py
List of Side Issues:
- [ ] Do not rely on
ValueErrorsin PyROS - [ ] Add custom callback to PyROS for debugging, etc.
Next steps from 2/18/2022 meeting:
- [x] Update pynumero_ASL, need new IDAES get extensions?
- Updated pynumero_ASL to version 2
- IDAES extensions build versions:
- Solvers: v2.4.4 20210524-2032 ubuntu2004
- Library: v2.4.4 20210524-2034 ubuntu2004
- [x] Post output shown in meeting today here with brief description on how it was generated
- [x] Run code at link provided in original post. Pyomo fork/branch: https://github.com/shermanjasonaf/pyomo/tree/debug-rankine-check-rank
- [x] Residuals of equations
- With tol of 1e-6, output message:
No constraints with residuals larger than 1e-06 !
- With tol of 1e-6, output message:
- [x] Errors when trying to check rank:
Traceback (most recent call last):
File "multiperiod_rankine_uncertainty.py", line 137, in <module>
results = pyros_solver.solve(
File "/home/jasherma/Documents/cmu/phd-project/pyomo_repo/pyomo/pyomo/contrib/pyros/pyros.py", line 421, in solve
pyros_soln, final_iter_separation_solns = ROSolver_iterative_solve(model_data, config)
File "/home/jasherma/Documents/cmu/phd-project/pyomo_repo/pyomo/pyomo/contrib/pyros/pyros_algorithm_methods.py", line 265, in ROSolver_iterative_solve
separation_problem_methods.solve_separation_problem(model_data=separation_data, config=config)
File "/home/jasherma/Documents/cmu/phd-project/pyomo_repo/pyomo/pyomo/contrib/pyros/separation_problem_methods.py", line 313, in solve_separation_problem
exit_separation_loop = solver_call_separation(model_data=model_data,
File "/home/jasherma/Documents/cmu/phd-project/pyomo_repo/pyomo/pyomo/contrib/pyros/separation_problem_methods.py", line 478, in solver_call_separation
dh.check_rank_equality_constraints()
File "/home/jasherma/Documents/cmu/phd-project/idaes-pse/idaes/core/util/model_diagnostics.py", line 221, in check_rank_equality_constraints
self.svd_analysis()
File "/home/jasherma/Documents/cmu/phd-project/idaes-pse/idaes/core/util/model_diagnostics.py", line 522, in svd_analysis
u, s, v = svds(self.jac_eq, k = n_sv, which='SM')
File "/home/jasherma/envs/rankine_env/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 1869, in svds
eigvals, eigvec = eigsh(XH_X, k=k, tol=tol ** 2, maxiter=maxiter,
File "/home/jasherma/envs/rankine_env/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 1690, in eigsh
params.iterate()
File "/home/jasherma/envs/rankine_env/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 570, in iterate
self._raise_no_convergence()
File "/home/jasherma/envs/rankine_env/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 376, in _raise_no_convergence
raise ArpackNoConvergence(msg % (num_iter, k_ok, self.k), ev, vec)
scipy.sparse.linalg.eigen.arpack.arpack.ArpackNoConvergence: ARPACK error -1: No convergence (2601 iterations, 0/10 eigenvectors converged)
@dallan-keylogic @Robbybp Do you have an enhancement for degeneracy hunter for an alternate SVD analysis? If so, can you share with @shermanjasonaf?
Notes from 3/4/2022:
- Output above is for after Ipopt terminates.
- Motivating question: Are any constraints infeasible at the initial point for the separation problem? If so, which constraints? How can we improve the initialization?
- Recall this problem as a special structure where the uncertainty only appears (linearly) in the objective. This means the constraints for the Rankine cycle should all be feasible at the initial point of the separation problem.
Next steps:
- [x] Post Ipopt output to here before changing the number of iterations. Post a link to the GitHub commit (hash) that corresponds to the code.
- [x] Set Ipopt iteraton to zero, post Ipopt output here
- [x] Post degeneracy hunter output for checking residuals with zero Ipopt iterations
- [ ] If there are any infeasible constraints at iteration 0, brainstorm how to improve initialization
- [ ] Try simplifying the decision rule. Attempt solving with PyROS decision rule zero. Post results here.
- [x] IPOPT output for PyROS subproblem as is; see:
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 855
Number of nonzeros in inequality constraint Jacobian.: 12
Number of nonzeros in Lagrangian Hessian.............: 189
Total number of variables............................: 354
variables with only lower bounds: 103
variables with lower and upper bounds: 6
variables with only upper bounds: 0
Total number of equality constraints.................: 348
Total number of inequality constraints...............: 14
inequality constraints with only lower bounds: 7
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 7
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.7500406e+02 1.00e-02 4.06e-12 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.7500406e+02 1.00e-04 1.01e-06 -1.0 1.00e-02 - 9.90e-01 9.90e-01h 1
2 1.7500406e+02 1.00e-06 1.69e-03 -1.0 1.00e-04 - 9.90e-01 9.90e-01h 1
3 1.7500406e+02 8.45e-07 9.05e+03 -1.0 1.34e-06 - 9.90e-01 8.99e-01h 1
4 1.7500406e+02 7.15e-07 5.34e+01 -1.0 3.50e-01 - 9.94e-01 1.00e+00f 1
5 1.7500406e+02 5.48e-06 1.94e+03 -2.5 3.70e-01 - 1.00e+00 1.22e-04h 14
6 1.7500406e+02 5.60e-01 2.83e-08 -2.5 9.10e-01 - 1.00e+00 1.00e+00s 22
7 1.7500406e+02 8.31e-02 1.30e+04 -5.7 1.04e+00 - 1.00e+00 8.61e-01h 1
8 1.7500406e+02 1.62e-04 3.10e+04 -5.7 3.12e-01 - 9.14e-01 1.00e+00h 1
9 1.7500406e+02 4.77e-07 6.88e-11 -5.7 1.62e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.7500406e+02 4.77e-07 1.84e+03 -8.6 5.80e-05 - 1.00e+00 1.22e-04h 14
11 1.7500406e+02 4.77e-07 1.83e+03 -8.6 3.21e+00 - 2.60e-03 2.60e-03s 14
12 1.7500406e+02 4.77e-07 1.82e+03 -8.6 1.90e-02 - 5.48e-03 5.48e-03s 14
13r 1.7500406e+02 4.77e-07 1.00e+03 -7.2 0.00e+00 - 0.00e+00 0.00e+00R 1
14r 1.7500406e+02 4.77e-07 9.91e+02 -7.2 3.35e-02 - 9.33e-03 9.90e-04f 1
15r 1.7500406e+02 5.96e-07 9.88e+02 -7.2 3.02e-02 - 3.24e-03 5.68e-03h 2
16r 1.7500406e+02 4.77e-07 9.73e+02 -7.2 3.21e-02 - 1.46e-02 8.35e-03h 2
17r 1.7500406e+02 4.77e-07 9.40e+02 -7.2 2.91e-02 - 3.39e-02 4.34e-02h 1
18r 1.7500406e+02 4.77e-07 8.01e+02 -7.2 2.20e-02 - 1.48e-01 1.39e-01H 1
19r 1.7500406e+02 4.77e-07 2.61e+02 -7.2 8.00e-03 - 7.41e-01 2.89e-06h 18
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 1.7500406e+02 4.77e-07 3.40e+02 -7.2 4.42e-03 - 7.23e-01 7.63e-06h 18
21r 1.7500406e+02 4.77e-07 2.29e+02 -7.2 3.43e-03 - 4.97e-01 7.63e-06h 18
22r 1.7500406e+02 2.98e-07 2.07e+02 -7.2 2.48e-03 - 3.17e-01 5.00e-01h 2
23r 1.7500406e+02 4.77e-07 2.04e+02 -7.2 1.65e-03 - 1.00e+00 2.50e-01h 3
24r 1.7500406e+02 3.58e-07 2.41e+02 -7.2 1.78e-03 - 1.00e+00 1.95e-03h 10
25r 1.7500406e+02 4.77e-07 1.62e+02 -7.2 1.76e-03 - 1.00e+00 2.50e-01h 3
26r 1.7500406e+02 4.77e-07 1.54e+02 -7.2 1.33e-03 - 1.00e+00 1.25e-01h 4
27r 1.7500406e+02 4.77e-07 1.47e+02 -7.2 1.15e-03 - 1.00e+00 1.95e-03h 10
28r 1.7500406e+02 4.77e-07 9.09e-13 -7.2 1.15e-03 - 1.00e+00 1.00e+00s 34
29r 1.7500406e+02 4.77e-07 1.19e+01 -9.0 2.75e-03 - 1.00e+00 8.57e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 1.7500406e+02 4.77e-07 2.25e+02 -9.0 2.67e-04 - 1.00e+00 2.74e-02f 6
31r 1.7500406e+02 4.77e-07 1.10e+01 -9.0 2.25e-04 - 9.51e-01 9.51e-01s 39
32r 1.7500406e+02 3.58e-07 2.27e-13 -9.0 4.06e-05 - 1.00e+00 1.00e+00s 39
Restoration phase converged to a point with small primal infeasibility.
Number of Iterations....: 32
(scaled) (unscaled)
Objective...............: 1.7500405522831917e+02 1.7500405522831917e+02
Dual infeasibility......: 2.9826141556554830e-13 2.9826141556554830e-13
Constraint violation....: 5.9604644775390625e-08 3.5762786865234375e-07
Complementarity.........: 9.0911772832484302e-10 9.0911772832484302e-10
Overall NLP error.......: 5.9604644775390625e-08 3.5762786865234375e-07
Number of objective function evaluations = 272
Number of objective gradient evaluations = 15
Number of equality constraint evaluations = 272
Number of inequality constraint evaluations = 272
Number of equality constraint Jacobian evaluations = 36
Number of inequality constraint Jacobian evaluations = 36
Number of Lagrangian Hessian evaluations = 33
Total CPU secs in IPOPT (w/o function evaluations) = 0.061
Total CPU secs in NLP function evaluations = 1.231
EXIT: Restoration Failed!
- [x] IPOPT
max_iter=0results (note: same traceback produced bycheck_rank_equality_constraints): - [x] PyROS separation subproblem solver Tee:
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 855
Number of nonzeros in inequality constraint Jacobian.: 12
Number of nonzeros in Lagrangian Hessian.............: 189
Total number of variables............................: 354
variables with only lower bounds: 103
variables with lower and upper bounds: 6
variables with only upper bounds: 0
Total number of equality constraints.................: 348
Total number of inequality constraints...............: 14
inequality constraints with only lower bounds: 7
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 7
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.7500406e+02 1.00e-02 4.06e-12 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
Number of Iterations....: 0
(scaled) (unscaled)
Objective...............: 1.7500405522831917e+02 1.7500405522831917e+02
Dual infeasibility......: 4.0575032881303200e-12 4.0575032881303200e-12
Constraint violation....: 9.9999916979527370e-03 9.9999916979527370e-03
Complementarity.........: 2.4230000000000011e+07 2.4230000000000011e+07
Overall NLP error.......: 2.4230000000000011e+07 2.4230000000000011e+07
Number of objective function evaluations = 1
Number of objective gradient evaluations = 1
Number of equality constraint evaluations = 1
Number of inequality constraint evaluations = 1
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 0
Total CPU secs in IPOPT (w/o function evaluations) = 0.020
Total CPU secs in NLP function evaluations = 0.000
EXIT: Maximum Number of Iterations Exceeded.
- [x] Degeneracy hunter check residuals output (for subproblem with
max_iter=0)
All constraints with residuals larger than 1e-05 :
count = 0 |residual| = 0.009999991697941343
energy_change : Size=1, Index=None, Active=True
Key : Lower : Body : Upper : Active
None : 0.0 : blocks[4].process.battery.next_soc - (blocks[4].process.battery.soc - blocks[4].process.battery.discharge + blocks[4].process.rankine.P_to_battery) : 0.0 : True
variable lower value upper
blocks[4].process.battery.next_soc 0 0.00999999 None
blocks[4].process.battery.soc 0 50.000000498363235 None
blocks[4].process.rankine.P_to_battery 0 224.99601701611587 None
blocks[4].process.battery.discharge None 274.99601751617706 None
- [ ] Enhanced decision rule order 0 simplification result in the works.
I was just dealing with a bunch of Restoration Failures with IAPWs 95. What @eslickj suggested is that they came from using saturated liquid or vapor enthalpies in constraints. See this constraint in the simple Rankine cycle file, for example:
m.fs.pre_condenser.eq_outlet_cond = Constraint(
expr=m.fs.pre_condenser.control_volume.properties_out[0].enth_mol
== m.fs.pre_condenser.control_volume.properties_out[
0
].enth_mol_sat_phase["Liq"]
)
The problem is that the behavior of water is nonsmooth at those points, and therefore the directional derivatives don't agree.
The way I'm dealing with it right now is by adding in a small margin of subcooling or superheating to those constraints, so that IPOPT doesn't land right on the phase transition. You want to do this directly with enthalpy, not with a property like temperature_saturation (because temperature's nonsmooth behavior wrt changes in enthalpy causes IPOPT a lot of problems trying to solve the KKT conditions):
m.fs.pre_condenser.eq_outlet_cond = Constraint(
expr=m.fs.pre_condenser.control_volume.properties_out[0].enth_mol
== 0.995*m.fs.pre_condenser.control_volume.properties_out[
0
].enth_mol_sat_phase["Liq"]
)
This isn't a great solution, and I'd like to figure out an easier solution for evaporators and condensers.
Also, while that might make the error message go away, it should be a warning sign that sensitivity analysis based on partial derivatives may not give you the whole story here. Some Monte Carlo simulations may be in order.
@shermanjasonaf At the initial point (iter=0), the most violated constraint is the energy balance on the battery. It looks like setting blocks[4].process.battery.next_soc to 0 would cause this constraint to be feasible. Some thoughts:
- One possibility is that the problem is correctly initialized, but the IPOPT bound push is adjusting the value of the variable
blocks[4].process.battery.next_socfrom0to0.00999999. To check this hypothesis, adjust thebound_pushsettings in IPOPT. I recommend trying1E-6and then rechecking the IPOPT output and maximum residual withmax_iter=0. - Another possibility is that the problem is not being initialized correctly. After trying the suggestion above, try setting the value of
blocks[4].process.battery.next_socto0and see what constraints are infeasible at iteration 0.
What constraint was most infeasible when IPOPT terminated with a restoration failure?
@dallan-keylogic Thank you for pointing this out. We should also loop in @jghouse88 who has extensive experience with this toy problem (simple Rankine cycle) we are using. My suggestion is to look at which constraint(s) are infeasible when IPOPT terminates and then go from there.
@adowling2 From the solver trace @shermanjasonaf posted, it looks like all constraints are solved to within 1e-6, so I don't think primal infeasibility is the driver here. An added complication is that Pyomo doesn't load solver results for a restoration failure, but you can work around that by setting iter to one iteration before the restoration failure happens.
A quick and dirty method to make IPOPT terminate is to choose which inequality constraints are active at a solution and turn them into equality constraints. I think it will complain about multiplier calculations failing, but for a square problem I'm pretty sure all it cares about is primal infeasibility.
@adowling2 The variable blocks[4].process.battery.next_soc is, in fact, initialized to 0 for the separation problem in question. With bound_push=1e-6, the DegeneracyHunter.check_residuals method reports no residuals larger than 1e-5 for the solution at iteration 0. The IPOPT output (for max_iter=0) is below.
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 855
Number of nonzeros in inequality constraint Jacobian.: 12
Number of nonzeros in Lagrangian Hessian.............: 189
Total number of variables............................: 354
variables with only lower bounds: 103
variables with lower and upper bounds: 6
variables with only upper bounds: 0
Total number of equality constraints.................: 348
Total number of inequality constraints...............: 14
inequality constraints with only lower bounds: 7
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 7
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.7500406e+02 9.92e-07 4.06e-12 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
Number of Iterations....: 0
(scaled) (unscaled)
Objective...............: 1.7500405522831917e+02 1.7500405522831917e+02
Dual infeasibility......: 4.0575032881303200e-12 4.0575032881303200e-12
Constraint violation....: 9.9169795930720284e-07 9.9169795930720284e-07
Complementarity.........: 2.4230000000000011e+07 2.4230000000000011e+07
Overall NLP error.......: 2.4230000000000011e+07 2.4230000000000011e+07
Number of objective function evaluations = 1
Number of objective gradient evaluations = 1
Number of equality constraint evaluations = 1
Number of inequality constraint evaluations = 1
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 0
Total CPU secs in IPOPT (w/o function evaluations) = 0.025
Total CPU secs in NLP function evaluations = 0.000
EXIT: Maximum Number of Iterations Exceeded.
In accordance with @dallan-keylogic 's suggestion, I have set max_iter=32 to capture the solution just before restoration failure; DegeneracyHunter.check_residuals reports that the most infeasible constraint is blocks[3].process.rankine.fs.turbine.isentropic_energy_balance[0.0], with a residual of 3.5762786865234375e-07.
I didn't read through everything, so I may have missed this, but have you considered scaling?
@eslickj Scaling of model components is currently not considered or applied by PyROS; I will look into scaling considerations.
@shermanjasonaf Are you sure about that? My understanding -- @eslickj and @andrewlee94 correct me please -- is that IDAES implements scaling by passing custom scaling to Ipopt. My first recommendation would be to post the options being passed to Ipopt when it is called by PyROS.
@adowling2 I mean to say that PyROS itself does not scale model components for any sub-problems; as I understand it, options passed to IPOPT are applied to PyROS subproblems at solve time (only when the IPOPT solver's solve method is invoked). IPOPT options currently used for the multiperiod rankine model implementation are the defaults. If I instead obtain an IPOPT solver through the method idaes.core.util.get_solver (with default arguments) the IPOPT options are nlp_scaling_method: 'gradient-based' and tol: 1e-06, and in attempting to solve with PyROS the same exceptions discussed above raised.
I'm super late, but IDAES scaling is slightly complicated. The way it works (very roughly) is that you specify scaling on a subset of variables. Usually like pressure, flow, sizes... maybe a few other things. IDAES tends to generate warnings about missing scaling factors if they are used but not supplied. From there you call calculate_scaling_factors(model) on the whole model. This uses the basic set of scaling factors to calculate scaling factors for constraints and intermediate variables.
So here's the tricky parts. The variable scale factors are ignored by IPOPT unless you set scaling_method to "user-scaling" (think I'm remembering the option right). Switching to user scaling causes whatever autoscaling ipopt does to be disabled. Depending on who you ask, you may or may not want the variable scaling to go to Ipopt and forgo the autoscaling. I do use user scaling.
Then for constraints things get extra complicated. Generally in IDAES we are transforming the constraints by just multiplying them by a scale factor. We usually consider that we want the constraint scaled so that a feasibility tolerance of like 1e-8 makes senses. There is another way to scale equations which is to supply a scaling factor. Ipopt's autoscale and the equivalent auto scaling in IDAES (rarely used) supply scaling factors for constraints. As I understand it, these scaling factors don't affect the tolerances used to check for feasibility, they just scale the Jacobian for the linear solver. The idea behind that is that whoever wrote the problem had the tolerances in mind when they wrote the constraints, and they wrote something that makes sense. That's not really true of IDAES, which is why we transform the equations based on the variable scaling so that tolerances make sense, and that is why you still may prefer to stick with ipopts gradient based scaling when you solve. So there are two different ways to scale constraints and they consider different things.
So all that is to say, as long as you call calculate_scale_factors(model), it still may help a lot. Whether you care about the scaling suffix is apparently debatable.
