Description
Summary
I’ve been experimenting with different solvers to see which ones to use under which robust optimizations. So, I started with the example provided under https://pyomo.readthedocs.io/en/stable/contributed_packages/pyros.html#pyros-usage-example. However, I was unable to get the single-stage version of this problem to work with IPOPT. The code in the example uses the baron
solver. However, since there don’t seem to be any integral uncertain parameters, nor integral variables, this problem should be solvable by a non-linear solver, e.g. IPOPT. I tried this example with IPOPT and got the following error:
Steps to reproduce the issue
# === Required import ===
import pyomo.environ as pyo
import pyomo.contrib.pyros as pyros
# === Instantiate the PyROS solver object ===
pyros_solver = pyo.SolverFactory("pyros")
# === Construct the Pyomo model object ===
m = pyo.ConcreteModel()
m.name = "hydro"
# === Define variables ===
m.x1 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
m.x2 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
m.x3 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
m.x4 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
m.x5 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
m.x6 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
m.x7 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
m.x8 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
m.x9 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
m.x10 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
m.x11 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
m.x12 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
m.x13 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x14 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x15 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x16 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x17 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x18 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x19 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x20 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x21 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x22 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x23 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x24 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
m.x25 = pyo.Var(within=pyo.Reals,bounds=(100000,100000),initialize=100000)
m.x26 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
m.x27 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
m.x28 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
m.x29 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
m.x30 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
m.x31 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
# === Define parameters ===
m.set_of_params = pyo.Set(initialize=[0, 1, 2, 3])
nominal_values = {0:82.8*0.0016, 1:4.97, 2:4.97, 3:1800}
m.p = pyo.Param(m.set_of_params, initialize=nominal_values, mutable=True)
# === Specify the objective function ===
m.obj = pyo.Objective(expr=m.p[0]*m.x1**2 + 82.8*8*m.x1 + 82.8*0.0016*m.x2**2 +
82.8*82.8*8*m.x2 + 82.8*0.0016*m.x3**2 + 82.8*8*m.x3 +
82.8*0.0016*m.x4**2 + 82.8*8*m.x4 + 82.8*0.0016*m.x5**2 +
82.8*8*m.x5 + 82.8*0.0016*m.x6**2 + 82.8*8*m.x6 + 248400,
sense=pyo.minimize)
# === Specify the constraints ===
m.c2 = pyo.Constraint(expr=-m.x1 - m.x7 + m.x13 + 1200<= 0)
m.c3 = pyo.Constraint(expr=-m.x2 - m.x8 + m.x14 + 1500 <= 0)
m.c4 = pyo.Constraint(expr=-m.x3 - m.x9 + m.x15 + 1100 <= 0)
m.c5 = pyo.Constraint(expr=-m.x4 - m.x10 + m.x16 + m.p[3] <= 0)
m.c6 = pyo.Constraint(expr=-m.x5 - m.x11 + m.x17 + 950 <= 0)
m.c7 = pyo.Constraint(expr=-m.x6 - m.x12 + m.x18 + 1300 <= 0)
m.c8 = pyo.Constraint(expr=12*m.x19 - m.x25 + m.x26 == 24000)
m.c9 = pyo.Constraint(expr=12*m.x20 - m.x26 + m.x27 == 24000)
m.c10 = pyo.Constraint(expr=12*m.x21 - m.x27 + m.x28 == 24000)
m.c11 = pyo.Constraint(expr=12*m.x22 - m.x28 + m.x29 == 24000)
m.c12 = pyo.Constraint(expr=12*m.x23 - m.x29 + m.x30 == 24000)
m.c13 = pyo.Constraint(expr=12*m.x24 - m.x30 + m.x31 == 24000)
m.c14 = pyo.Constraint(expr=-8e-5*m.x7**2 + m.x13 == 0)
m.c15 = pyo.Constraint(expr=-8e-5*m.x8**2 + m.x14 == 0)
m.c16 = pyo.Constraint(expr=-8e-5*m.x9**2 + m.x15 == 0)
m.c17 = pyo.Constraint(expr=-8e-5*m.x10**2 + m.x16 == 0)
m.c18 = pyo.Constraint(expr=-8e-5*m.x11**2 + m.x17 == 0)
m.c19 = pyo.Constraint(expr=-8e-5*m.x12**2 + m.x18 == 0)
m.c20 = pyo.Constraint(expr=-4.97*m.x7 + m.x19 == 330)
m.c21 = pyo.Constraint(expr=-m.p[1]*m.x8 + m.x20 == 330)
m.c22 = pyo.Constraint(expr=-4.97*m.x9 + m.x21 == 330)
m.c23 = pyo.Constraint(expr=-4.97*m.x10 + m.x22 == 330)
m.c24 = pyo.Constraint(expr=-m.p[2]*m.x11 + m.x23 == 330)
m.c25 = pyo.Constraint(expr=-4.97*m.x12 + m.x24 == 330)
# === Specify which parameters are uncertain ===
uncertain_parameters = [m.p] # We can pass IndexedParams this way to PyROS, or as an expanded list per index
# === Define the pertinent data ===
relative_deviation = 0.15
bounds = [(nominal_values[i] - relative_deviation*nominal_values[i],
nominal_values[i] + relative_deviation*nominal_values[i])
for i in range(4)]
# === Construct the desirable uncertainty set ===
box_uncertainty_set = pyros.BoxSet(bounds=bounds)
ipopt_solver = pyo.SolverFactory('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.x1, m.x2, m.x3, m.x4, m.x5, m.x6,
m.x19, m.x20, m.x21, m.x22, m.x23, m.x24, m.x31]
second_stage_variables = []
# The remaining variables are implicitly designated to be state variables
# === Call PyROS to solve the robust optimization problem ===
results_1 = pyros_solver.solve(model = m,
first_stage_variables = first_stage_variables,
second_stage_variables = second_stage_variables,
uncertain_params = uncertain_parameters,
uncertainty_set = box_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
})
# === Query results ===
time = results_1.time
iterations = results_1.iterations
termination_condition = results_1.pyros_termination_condition
objective = results_1.final_objective_value
# === Print some results ===
single_stage_final_objective = round(objective,-1)
print("Final objective value: %s" % single_stage_final_objective)
print("PyROS termination condition: %s" % termination_condition)
Error Message
Here's the error message:
===========================================================================================
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.
===========================================================================================
WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
- termination condition: other
- message from solver: <undefined>
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Input In [5], in <cell line: 107>()
103 second_stage_variables = []
104 # The remaining variables are implicitly designated to be state variables
105
106 # === Call PyROS to solve the robust optimization problem ===
--> 107 results_1 = pyros_solver.solve(model = m,
108 first_stage_variables = first_stage_variables,
109 second_stage_variables = second_stage_variables,
110 uncertain_params = uncertain_parameters,
111 uncertainty_set = box_uncertainty_set,
112 local_solver = local_solver,
113 global_solver= global_solver,
114 tee=True,
115 options = {
116 "objective_focus": pyros.ObjectiveType.worst_case,
117 "solve_master_globally": True,
118 "load_solution":False
119 })
121 # === Query results ===
122 time = results_1.time
File ~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/pyros.py:365, in PyROS.solve(self, model, first_stage_variables, second_stage_variables, uncertain_params, uncertainty_set, local_solver, global_solver, **kwds)
361 model.add_component(model_data.util_block, util)
362 # Note: model.component(model_data.util_block) is util
363
364 # === Validate uncertainty set happens here, requires util block for Cardinality and FactorModel sets
--> 365 validate_uncertainty_set(config=config)
367 # === Leads to a logger warning here for inactive obj when cloning
368 model_data.original_model = model
File ~/opt/anaconda3/lib/python3.9/site-packages/pyomo/contrib/pyros/util.py:211, in validate_uncertainty_set(config)
209 # === Check set validity via boundedness and non-emptiness
210 if not config.uncertainty_set.is_valid(config=config):
--> 211 raise AttributeError("Invalid uncertainty set detected. Check the uncertainty set object to "
212 "ensure non-emptiness and boundedness.")
214 return
AttributeError: Invalid uncertainty set detected. Check the uncertainty set object to ensure non-emptiness and boundedness.
Information on your system
Pyomo version: 6.4.1
Python version: 3.9
Operating system: MacOS Monterey
Solver (if applicable): IPOPT
Additional information
Is this the expected behavior? Or am I failing to see how the problem exceeds the capabilities of a non-linear solver such as IPOPT?