Pyomo
Nominal value affects robust feasibility even when `solve_master_globally=True`?
Summary
Below I’ve shown a very simple problem for which PyROS is still unable to find a robust feasible solution when I use IPOPT as the local and global solver. However, if the nominal value is changed, then PyROS is able to find a robust feasible solution.
Notice that z=0 is feasible for any realization of the uncertain parameters. This might explain my difficulty in getting PyROS to return a robust feasible solution when a similar problem is solved on a larger dimension. What effect does the nominal value have on robust feasibility?
The only information I saw on this page https://pyomo.readthedocs.io/en/stable/contributed_packages/pyros.html was the following:
“If objective_focus = ObjectiveType.nominal, second-stage objective and variables are evaluated at the nominal realization of the uncertain parameters...” However, for several reasons, this information doesn’t apply to the problem below.
As a side note, if the objective function is changed such that it involves the uncertain parameter “d”, then PyROS seems to find a robust feasible solution regardless of its (feasible) nominal value.
And how should I set my nominal values to help PyROS find this robust feasible solution?
Steps to reproduce the issue
import numpy as np
def create_model():
m = pe.ConcreteModel()
m.del_component('q')
m.add_component('q', pe.Param(mutable=True, initialize=0.5))
m.del_component('d')
# If we change the initial value to 1.0 instead of 0.0, a robust feasible solution is found
m.add_component('d', pe.Param(mutable=True, initialize=0.0))
m.del_component('z')
m.add_component('z', pe.Var(within=pe.Reals, initialize=1))
return m
m = create_model()
m.constraints_expression_1 = pe.Constraint(expr = 0<= m.z*m.d*m.q)
# === Specify the objective function ===
m.del_component('objective_fn_expression')
m.add_component('objective_fn_expression', pe.Objective(expr = m.z*m.q, sense = pe.minimize))
import pyomo.contrib.pyros as pyros
# === Instantiate the PyROS solver object ===
pyros_solver = pe.SolverFactory("pyros")
# === Specify which parameters are uncertain ===
uncertain_parameters = [m.d, m.q]
lhs_coefficients_mat = np.array([[1.0,0.0],[-1.0,0.0], [0.0,1.0],[0.0,-1.0]])
rhs_vec = np.array([1.0,0.0,1000.0,0.0])
uncertainty_set = pyros.PolyhedralSet(lhs_coefficients_mat, rhs_vec)
# === Designate local and global NLP solvers ===
ipopt_solver = pe.SolverFactory('appsi_ipopt')
local_solver = ipopt_solver
global_solver = ipopt_solver
# === Designate which variables correspond to first- and second-stage degrees of freedom ===
first_stage_variables = [m.z]
second_stage_variables = []
# The remaining variables are implicitly designated to be state variables
# === Call PyROS to solve the robust optimization problem ===
results = pyros_solver.solve(model = m,
first_stage_variables = first_stage_variables,
second_stage_variables = second_stage_variables,
uncertain_params = uncertain_parameters,
uncertainty_set = uncertainty_set,
local_solver = local_solver,
global_solver = global_solver,
tee = False,
options = {
"objective_focus": pyros.ObjectiveType.worst_case,
"solve_master_globally": True,
"load_solution":False,
"bypass_global_separation":False
})
# === Query results ===
time = results.time
iterations = results.iterations
termination_condition = results.pyros_termination_condition
objective = results.final_objective_value
# === Print some results ===
RuntimeError: A feasible solution was not found, so no solution can be loaded.Please set opt.config.load_solution=False and check results.termination_condition and resutls.best_feasible_objective before loading a solution.
@makansij Your deterministic model (and therefore the initial master problem) is unbounded when d is 0 and q is nonzero, since there are no explicit variable bounds for z and your objective expression is z * q. PyROS doesn't properly account for the exception raised by your subsolver in this case. (Run PyROS with tee=True to see the subsolver output.)
The nominal uncertain parameter realization should not affect the robust feasibility of your model, but it may indeed affect the trajectory of the iterates, and therefore possibly the termination outcome. Moreover, the final PyROS outcome may change if there is an issue solving a subproblem successfully (such as a local infeasibility, or some other numerical/convergence issue).
I would verify that at least the local subsolver is capable of solving your determinstic model to optimality (under at least your nominal parameter realization) before using PyROS.
@shermanjasonaf , below is an example of a problem for which the deterministic master problem can be solved at nominal values, but the problem with uncertain parameters cannot. I can think of at least one robust feasible solution (the trivial one where Z is all 0's).
You'll notice, if you flip the uncertain_flag, then the problem is solved.
import pyomo.environ as pe
import scipy.io as sp_io
import pyomo.contrib.pyros as pyros
import numpy as np
def setup_model():
m = pe.ConcreteModel()
m.del_component('N')
m.add_component('N', pe.Param(mutable=True, within=pe.Reals, initialize=2))
N = pe.value(m.N)
m.del_component('n')
m.add_component('n', pe.Param(mutable=True, within=pe.Reals, initialize=3))
n = pe.value(m.n)
def y_init(m, i):
y_val = 0.5*np.ones((1,pe.value(m.n)))
return y_val[0,i]
m.del_component('y')
m.add_component('y', pe.Param(range(n), mutable=True, within=pe.Reals, initialize=y_init))
m.del_component('Q')
m.add_component('Q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=1))
Q = pe.value(m.Q)
m.del_component('q')
m.add_component('q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=0))
def d_init(m, ix):
return float(1/m.Q.value)
m.del_component('d')
m.add_component('d', pe.Param(range(m.Q.value), mutable=True, within=pe.Reals, initialize=d_init))
def r_init(m, i):
r = np.array([[1,0,0]])
return r[0, i]
m.del_component('r')
m.add_component('r', pe.Param(range(n), initialize=r_init)) # 1 by n
# Finally, add the decision variable
m.del_component('Z')
m.add_component('Z', pe.Var(range(n), range(n*Q*N), within=pe.Reals))
return m
m = setup_model()
m.del_component('constraints_expression_1_0_0')
m.add_component('constraints_expression_1_0_0', pe.Constraint(expr = 0<= m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2]))
m.del_component('constraints_expression_1_0_1')
m.add_component('constraints_expression_1_0_1', pe.Constraint(expr = 0<= m.y[0]*m.d[0]*m.Z[0,3] + m.y[1]*m.d[0]*m.Z[0,4] + m.d[0]*m.Z[0,5]*m.y[2]))
m.del_component('constraints_expression_1_1_0')
m.add_component('constraints_expression_1_1_0', pe.Constraint(expr = 0<= m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2]))
m.del_component('constraints_expression_1_1_1')
m.add_component('constraints_expression_1_1_1', pe.Constraint(expr = 0<= m.y[0]*m.d[0]*m.Z[1,3] + m.y[1]*m.d[0]*m.Z[1,4] + m.d[0]*m.Z[1,5]*m.y[2]))
m.del_component('constraints_expression_1_2_0')
m.add_component('constraints_expression_1_2_0', pe.Constraint(expr = 0<= m.y[0]*m.d[0]*m.Z[2,0] + m.y[1]*m.d[0]*m.Z[2,1] + m.d[0]*m.Z[2,2]*m.y[2]))
m.del_component('constraints_expression_1_2_1')
m.add_component('constraints_expression_1_2_1', pe.Constraint(expr = 0<= m.y[0]*m.d[0]*m.Z[2,3] + m.y[1]*m.d[0]*m.Z[2,4] + m.d[0]*m.Z[2,5]*m.y[2]))
m.del_component('constraints_expression_2_0_0')
m.add_component('constraints_expression_2_0_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2] <= 0.9*m.y[0]))
m.del_component('constraints_expression_2_0_1')
m.add_component('constraints_expression_2_0_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,3] + m.y[1]*m.d[0]*m.Z[0,4] + m.d[0]*m.Z[0,5]*m.y[2] - (0.9*m.y[0] + 0.9 - 0.9*(m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.y[2]*m.d[0]*m.Z[0,2] )) <=0 ))
m.del_component('constraints_expression_2_1_0')
m.add_component('constraints_expression_2_1_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2] <= 0.9*m.y[1]))
m.del_component('constraints_expression_2_1_1')
m.add_component('constraints_expression_2_1_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,3] + m.y[1]*m.d[0]*m.Z[1,4] + m.d[0]*m.Z[1,5]*m.y[2] - (0.9*m.y[1] + (m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2]) - 0.9*(m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.y[2]*m.d[0]*m.Z[1,2])) <= 0))
m.del_component('constraints_expression_2_2_0')
m.add_component('constraints_expression_2_2_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[2,0] + m.y[1]*m.d[0]*m.Z[2,1] + m.d[0]*m.Z[2,2]*m.y[2] <= 0.9*m.y[2]))
m.del_component('constraints_expression_2_2_1')
m.add_component('constraints_expression_2_2_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[2,3] + m.y[1]*m.d[0]*m.Z[2,4] + m.d[0]*m.Z[2,5]*m.y[2] - (0.9*m.y[2] + (m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2]) - 0.9*(m.y[0]*m.d[0]*m.Z[2,0] + m.y[1]*m.d[0]*m.Z[2,1] + m.y[2]*m.d[0]*m.Z[2,2])) <= 0.0))
m.del_component('constraints_expression_3_0_0')
m.add_component('constraints_expression_3_0_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2] <= 100.0))
m.del_component('constraints_expression_3_0_1')
m.add_component('constraints_expression_3_0_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,3] + m.y[1]*m.d[0]*m.Z[0,4] + m.d[0]*m.Z[0,5]*m.y[2] <= 100.0 ))
m.del_component('constraints_expression_3_1_0')
m.add_component('constraints_expression_3_1_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2] <= 100.0))
m.del_component('constraints_expression_3_1_1')
m.add_component('constraints_expression_3_1_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,3] + m.y[1]*m.d[0]*m.Z[1,4] + m.d[0]*m.Z[1,5]*m.y[2]<= 100.0))
m.del_component('constraints_expression_3_2_0')
m.add_component('constraints_expression_3_2_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[2,0] + m.y[1]*m.d[0]*m.Z[2,1] + m.d[0]*m.Z[2,2]*m.y[2]<= 100.0))
m.del_component('constraints_expression_3_2_1')
m.add_component('constraints_expression_3_2_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[2,3] + m.y[1]*m.d[0]*m.Z[2,4] + m.d[0]*m.Z[2,5]*m.y[2]<= 100.0))
m.del_component('constraints_expression_4_0_1_0')
m.add_component('constraints_expression_4_0_1_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2] <= 1 - 0.3*m.y[1]))
m.del_component('constraints_expression_4_0_1_1')
m.add_component('constraints_expression_4_0_1_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,3] + m.y[1]*m.d[0]*m.Z[0,4] + m.d[0]*m.Z[0,5]*m.y[2] - (1 - (0.3*m.y[1] + 0.3*(m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2]) - 0.3*(m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.y[2]*m.d[0]*m.Z[1,2]))) <= 0.0))
m.del_component('constraints_expression_4_1_2_0')
m.add_component('constraints_expression_4_1_2_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2] <= 1 - 0.3*m.y[2]))
m.del_component('constraints_expression_4_1_2_1')
m.add_component('constraints_expression_4_1_2_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,3] + m.y[1]*m.d[0]*m.Z[1,4] + m.d[0]*m.Z[1,5]*m.y[2] - (1 - (0.3*m.y[2] + 0.3*(m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2]) - 0.3*(m.y[0]*m.d[0]*m.Z[2,0] + m.y[1]*m.d[0]*m.Z[2,1] + m.y[2]*m.d[0]*m.Z[2,2]))) <= 0.0))
m.del_component('constraints_expression_5_0_1_0')
m.add_component('constraints_expression_5_0_1_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,0] + m.y[1]*m.d[0]*m.Z[0,1] + m.d[0]*m.Z[0,2]*m.y[2] <= 1))
m.del_component('constraints_expression_5_0_1_1')
m.add_component('constraints_expression_5_0_1_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[0,3] + m.y[1]*m.d[0]*m.Z[0,4] + m.d[0]*m.Z[0,5]*m.y[2] <= 1))
m.del_component('constraints_expression_5_1_2_0')
m.add_component('constraints_expression_5_1_2_0', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,0] + m.y[1]*m.d[0]*m.Z[1,1] + m.d[0]*m.Z[1,2]*m.y[2] <= 1))
m.del_component('constraints_expression_5_1_2_1')
m.add_component('constraints_expression_5_1_2_1', pe.Constraint(expr = m.y[0]*m.d[0]*m.Z[1,3] + m.y[1]*m.d[0]*m.Z[1,4] + m.d[0]*m.Z[1,5]*m.y[2] <= 1))
uncertain_flag = True
# === Specify the objective function ===
m.del_component('objective_fn_expression')
m.add_component('objective_fn_expression', pe.Objective(expr = 4.5 - (4.5*m.Z[0,0] + 1.5*m.Z[0,3] - 3.0*m.Z[1,0] - 2.0*m.Z[1,3] - 3.0*m.Z[2,0] - m.Z[2,3] + (9*m.Z[0,1] + 3*m.Z[0,4] - 6*m.Z[1,1] - 4*m.Z[1,4] - 6*m.Z[2,1] - 2*m.Z[2,4])*0.5 + (9*m.Z[0,2] + 3*m.Z[0,5] - 6*m.Z[1,2] - 4*m.Z[1,5] - 6*m.Z[2,2] - 2*m.Z[2,5])*0.5), sense = pe.minimize))
if uncertain_flag:
import pyomo.contrib.pyros as pyros
# === Instantiate the PyROS solver object ===
pyros_solver = pe.SolverFactory("pyros")
# === Specify which parameters are uncertain ===
uncertain_parameters = [m.d, m.y]
lhs_coefficients_mat = np.array([[ 1., 0., 0., 0.], [-1., 0., 0., 0.], [ 1., 0., 0., 0.], [-1., 0., 0., 0.], [ 0., 1., 0., 0.], [ 0., -1., 0., 0.], [ 0., 0., 1., 0.], [ 0., 0., -1., 0.], [ 0., 0., 0., 1.], [ 0., 0., 0., -1.]])
rhs_vec = np.array([1.0001, -0.9999, 1.0, 0.0, 10.0, 0.0, 3.3333333333333335, 0.0, 3.3333333333333335, 0.0])
uncertainty_set = pyros.PolyhedralSet(lhs_coefficients_mat, rhs_vec)
# === Designate local and global NLP solvers ===
ipopt_solver = pe.SolverFactory('appsi_ipopt')
local_solver = ipopt_solver
global_solver = ipopt_solver
# === Designate which variables correspond to first- and second-stage degrees of freedom ===
first_stage_variables = [m.Z]
second_stage_variables = []
# The remaining variables are implicitly designated to be state variables
# === Call PyROS to solve the robust optimization problem ===
results = pyros_solver.solve(model = m,
first_stage_variables = first_stage_variables,
second_stage_variables = second_stage_variables,
uncertain_params = uncertain_parameters,
uncertainty_set = uncertainty_set,
local_solver = local_solver,
global_solver = global_solver,
tee = True,
options = {
"objective_focus": pyros.ObjectiveType.worst_case,
"solve_master_globally": True,
"load_solution":False,
"bypass_global_separation":False
})
else:
solver = pe.SolverFactory('appsi_ipopt')
results = solver.solve(m)
This example might give insight into whether the solver is the culprit here, or something else that I'm missing about the model setup. Is there something I'm missing about how to get PyROS to find this robust feasible solution of m.Z being all 0's? Or should I try COUENNE?
Thanks again for your continued assistance on this.
@makansij Can you share the traceback you encounter when you attempt to solve with uncertain_flag=True (i.e. attempting to solve with PyROS)?
Before we continue, I want to question what I said early about “the deterministic master problem can be solved at nominal values”. The output from the deterministic master problem, as you'll discover is the following:
{'Problem': [{'Lower bound': -inf, 'Upper bound': -1.6950002669605402, 'Number of objectives': 1, 'Number of constraints': 0, 'Number of variables': 0, 'Sense': 1}], 'Solver': [{'Status': 'ok', 'Termination condition': 'optimal', 'Termination message': 'TerminationCondition.optimal'}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
Although it returns ‘Termination condition’: ‘optimal’, it also returns 'Lower bound': -inf. This is beyond my knowledge of IPOPT/Pyomo. The optimal objective value is indeed -1.695... but it’s not clear to me whether IPOPT was able to find this.
For what it’s worth, solving the deterministic master problem with glpk also gives ‘Termination condition’: ‘optimal’, but instead of 'Lower bound': -inf, the Lower bound and 'Upper bound' are the same correct objective function value of -1.695.... Am I understanding the output from IPOPT/PYomo correctly, and that I should ignore the fact that 'Lower bound': -inf?
Your deterministic problem is an LP with a minimization objective. An LP solver such as glpk can certify that the optimal solution is, in fact, globally optimal; correspondingly, the 'Lower bound' and 'Upper bound' fields are set to the value of your objective function expression (-1.695...). The local NLP solver 'appsi_ipopt' provides no global optimality guarantees, (though theoretically, the locally optimal solution found should be globally optimal since your model is linear---and therefore convex,) and provides no other information about the objective lower bound. For this reason, the 'Lower bound' field has been set to -inf.
Sure. Here's the tail end of the output from the solver iterations, and the traceback:
Number of Iterations....: 65
(scaled) (unscaled)
Objective...............: -1.6589180375602830e+00 -1.6589180375602830e+00
Dual infeasibility......: 9.9999929814892730e-01 9.9999929814892730e-01
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 6.4333355613283178e-09 6.4333355613283178e-09
Overall NLP error.......: 9.9999929814892730e-01 9.9999929814892730e-01
Number of objective function evaluations = 259
Number of objective gradient evaluations = 27
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 260
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 70
Number of Lagrangian Hessian evaluations = 66
Total seconds in IPOPT = 0.092
EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/var/folders/1g/pjkrsvtx455c4cm1_j7xrnz00000gp/T/ipykernel_10508/3883132644.py in <module>
24
25 # === Call PyROS to solve the robust optimization problem ===
---> 26 results = pyros_solver.solve(model = m,
27 first_stage_variables = first_stage_variables,
28 second_stage_variables = second_stage_variables,
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/pyros.py in solve(self, model, first_stage_variables, second_stage_variables, uncertain_params, uncertainty_set, local_solver, global_solver, **kwds)
431
432 # === Solve and load solution into model
--> 433 pyros_soln, final_iter_separation_solns = ROSolver_iterative_solve(model_data, config)
434
435
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/pyros_algorithm_methods.py in ROSolver_iterative_solve(model_data, config)
161 # === Solve Master Problem
162 config.progress_logger.info("PyROS working on iteration %s..." % k)
--> 163 master_soln = master_problem_methods.solve_master(model_data=master_data, config=config)
164 #config.progress_logger.info("Done solving Master Problem!")
165 master_soln.master_problem_subsolver_statuses = []
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/master_problem_methods.py in solve_master(model_data, config)
522 solver = config.global_solver if config.solve_master_globally else config.local_solver
523
--> 524 return solver_call_master(model_data=model_data, config=config, solver=solver,
525 solve_data=master_soln)
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/master_problem_methods.py in solver_call_master(model_data, config, solver, solve_data)
466 solver = backup_solvers.pop(0)
467 try:
--> 468 results = solver.solve(nlp_model, tee=config.tee)
469 except ValueError as err:
470 if 'Cannot load a SolverResults object with bad status: error' in str(err):
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/base.py in solve(self, model, tee, load_solutions, logfile, solnfile, timelimit, report_timing, solver_io, suffixes, options, keepfiles, symbolic_solver_labels)
1234 self.options = options
1235
-> 1236 results: Results = super(LegacySolverInterface, self).solve(model)
1237
1238 legacy_results = LegacySolverResults()
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/solvers/ipopt.py in solve(self, model, timer)
264 self._writer.write(model, self._filename+'.nl', timer=timer)
265 timer.stop('write nl file')
--> 266 res = self._apply_solver(timer)
267 self._last_results_object = res
268 if self.config.report_timing:
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/solvers/ipopt.py in _apply_solver(self, timer)
432 else:
433 timer.start('parse solution')
--> 434 results = self._parse_sol()
435 timer.stop('parse solution')
436
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/solvers/ipopt.py in _parse_sol(self)
371 results.best_feasible_objective = value(obj_expr_evaluated)
372 elif self.config.load_solution:
--> 373 raise RuntimeError('A feasible solution was not found, so no solution can be loaded.'
374 'Please set opt.config.load_solution=False and check '
375 'results.termination_condition and '
RuntimeError: A feasible solution was not found, so no solution can be loaded.Please set opt.config.load_solution=False and check results.termination_condition and resutls.best_feasible_objective before loading a solution
Re: glpk, thanks for your response. Got it.
Are you able to see that the value (the trivial one where Z is all 0's) for m.Z is robustly feasible?
I've run your script independently on my machine; PyROS converges to a robust feasible solution (in which m.Z is effectively all zeros) in five iterations. What versions of Pyomo and IPOPT are you using? Also, in which PyROS iteration do you encounter this traceback (run with tee=False)?
Interesting! Versions:
Ipopt 3.14.9, running with linear solver MUMPS 5.5.0. pyomo: 6.4.1
Do you mind sharing your versions ?
When tee=False, it fails at iteration 1. This is the output:
===========================================================================================
PyROS: Pyomo Robust Optimization Solver v.1.1.1
Developed by Natalie M. Isenberg (1), John D. Siirola (2), Chrysanthos E. Gounaris (1)
(1) Carnegie Mellon University, Department of Chemical Engineering
(2) Sandia National Laboratories, Center for Computing Research
The developers gratefully acknowledge support from the U.S. Department of Energy's
Institute for the Design of Advanced Energy Systems (IDAES)
===========================================================================================
======================================== DISCLAIMER =======================================
PyROS is still under development.
Please provide feedback and/or report any issues by opening a Pyomo ticket.
===========================================================================================
INFO: PyROS working on iteration 0...
INFO: PyROS working on iteration 1...
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/var/folders/1g/pjkrsvtx455c4cm1_j7xrnz00000gp/T/ipykernel_12283/517130743.py in <module>
24
25 # === Call PyROS to solve the robust optimization problem ===
---> 26 results = pyros_solver.solve(model = m,
27 first_stage_variables = first_stage_variables,
28 second_stage_variables = second_stage_variables,
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/pyros.py in solve(self, model, first_stage_variables, second_stage_variables, uncertain_params, uncertainty_set, local_solver, global_solver, **kwds)
431
432 # === Solve and load solution into model
--> 433 pyros_soln, final_iter_separation_solns = ROSolver_iterative_solve(model_data, config)
434
435
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/pyros_algorithm_methods.py in ROSolver_iterative_solve(model_data, config)
161 # === Solve Master Problem
162 config.progress_logger.info("PyROS working on iteration %s..." % k)
--> 163 master_soln = master_problem_methods.solve_master(model_data=master_data, config=config)
164 #config.progress_logger.info("Done solving Master Problem!")
165 master_soln.master_problem_subsolver_statuses = []
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/master_problem_methods.py in solve_master(model_data, config)
522 solver = config.global_solver if config.solve_master_globally else config.local_solver
523
--> 524 return solver_call_master(model_data=model_data, config=config, solver=solver,
525 solve_data=master_soln)
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/master_problem_methods.py in solver_call_master(model_data, config, solver, solve_data)
466 solver = backup_solvers.pop(0)
467 try:
--> 468 results = solver.solve(nlp_model, tee=config.tee)
469 except ValueError as err:
470 if 'Cannot load a SolverResults object with bad status: error' in str(err):
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/base.py in solve(self, model, tee, load_solutions, logfile, solnfile, timelimit, report_timing, solver_io, suffixes, options, keepfiles, symbolic_solver_labels)
1234 self.options = options
1235
-> 1236 results: Results = super(LegacySolverInterface, self).solve(model)
1237
1238 legacy_results = LegacySolverResults()
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/solvers/ipopt.py in solve(self, model, timer)
264 self._writer.write(model, self._filename+'.nl', timer=timer)
265 timer.stop('write nl file')
--> 266 res = self._apply_solver(timer)
267 self._last_results_object = res
268 if self.config.report_timing:
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/solvers/ipopt.py in _apply_solver(self, timer)
432 else:
433 timer.start('parse solution')
--> 434 results = self._parse_sol()
435 timer.stop('parse solution')
436
~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/appsi/solvers/ipopt.py in _parse_sol(self)
371 results.best_feasible_objective = value(obj_expr_evaluated)
372 elif self.config.load_solution:
--> 373 raise RuntimeError('A feasible solution was not found, so no solution can be loaded.'
374 'Please set opt.config.load_solution=False and check '
375 'results.termination_condition and '
RuntimeError: A feasible solution was not found, so no solution can be loaded.Please set opt.config.load_solution=False and check results.termination_condition and resutls.best_feasible_objective before loading a solution.
IPOPT 3.13.2, with linear solver MA27. Pyomo 6.4.1
I'll try those versions - thanks!
@shermanjasonaf Do you mind trying my versions to see if you get the same result?
I’m having some trouble getting IPOPT to read the dynamic library files from MA27 (part of HSL). Since MA27 is the default option for HSL, I tried just downloading HSL and linking it to IPOPT.
I first decided to try getting this working for the deterministic model (uncertain_flag=False), so I’m not even using PyROS in the following steps.
Here’s what I did:
I followed the directions here https://coin-or.github.io/Ipopt/INSTALL.html under HSL header
The output from the final command sudo make install indicates that /usr/local/lib/libcoinhsl.2.dylib is installed.
However, when I run the above code, changing the solver from solver = pe.SolverFactory(glpk) to solver = pe.SolverFactory('ipopt'), IPOPT is unable to find the dylib file:
Ipopt 3.14.9:
Exception of type: DYNAMIC_LIBRARY_FAILURE in file "Common/IpLibraryLoader.cpp" at line 67:
Exception message: dlopen(libhsl.dylib, 0x0002): tried: 'libhsl.dylib' (no such file), '/usr/local/lib/libhsl.dylib' (no such file),...
After following the installation steps, I didn’t see any files named libhsl.dylib but there was a file named libcoinhsl.dylib inside Third-Party-HSL. Upon reading this https://discourse.julialang.org/t/needs-help-on-ipopt-coinhsl-installation/40583/2 I tried to rename the file libcoinhsl.dylib inside Third-Party-HSL to libhsl.dylib. Other sources (https://stackoverflow.com/questions/58305144/trying-to-compile-hsl-to-get-ipopt ) also suggest manual renaming, which I’m always highly skeptical of. However, this causes other issues
Exception of type: DYNAMIC_LIBRARY_FAILURE in file "Common/IpLibraryLoader.cpp" at line 175:
Exception message: dlsym(0x215515860, MA27AD): symbol not found
I’m not able to get much information on this error.
According to the HSL docs, I shouldn’t need to recompile IPOPT in order to link it with HSL. So, I also tried setting the options such that IPOPT picks up the dynamic library:
# solver.options['OF_hsllib'] = 'libhsl.dylib'
# solver.options['OF_linear_solver'] = 'ma27'
But the same error persists. Is there anything you might suggest for running the same solver as you and generating the solution using HSL/MA27?
IPOPT 3.13.2, with linear solver MA27. Pyomo 6.4.1
I was finally able to get MA27 to work with IPOPT version 3.14.9, and I can verify that a Robust Feasible solution was found. It seems like the solver makes all the difference.
@makansij We find this issue particularly interesting: your original script exemplifies the case in which your determinstic model (master 0) is unbounded for the parameter realizations in some subset of your uncertainty set. We are thinking about best ways for PyROS to recover from such cases. In the meantime, proceed with using a different nominal point as discussed previously.