pyomo
pyomo copied to clipboard
Pyomo
How to help PyROS verify nonempty Intersection Set, when nominal point is clearly in the intersection?
I’m experimenting with intersection uncertainty sets, and I’ve run into a problem where:
-
The nominal point is in the intersection set
-
But PyROS implies that the intersection is empty when solving the
nlpin line 1094 ofuncertainty_sets.pyin the pyros library.
I’ve tried using a variety of solvers, and I’ve verified that point_in_set returns True when executed on the nominal point. Also, my problem is a general NLP, so the solvers that I’ve tried should be able to solve (IPOPT, COUENNE). The code is shown below.
There is an error from the solver when I run with tee=True, however it doesn’t tell me much, and the solver runs successfully on other problems.
ERROR: Solver log: Ipopt 3.14.9: option_file_name=/var/folders/1g/pjkrsvtx455c
4cm1_j7xrnz00000gp/T/tmpzni2ajgz_ipopt.opt bad line 121 of
/var/folders/1g/pjkrsvtx455c4cm1_j7xrnz00000gp/T/tmpkdxumge4.pyomo.nl: 1
array([-10000.]) libc++abi: terminating with uncaught exception of type
Ipopt::TNLP::INVALID_TNLP
The minimal example shown below has a flag MA_27_flag to be able to toggle between different solvers.
Is there another way to set up the intersection set that can get around this problem?
Thanks.
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 x_0_init(m, i):
y = 0.5*np.ones((1,3))
return y[0,i]
m.del_component('y')
m.add_component('y', pe.Param(range(n), mutable=True, initialize=x_0_init))
def alpha_init(m, k, i):
alpha = np.array([[0, 0, 0], [2,1,0], [3,2,1]])
return alpha[k, i]
m.del_component('alpha')
m.add_component('alpha', pe.Param(range(N+1), range(n), mutable=True, initialize=alpha_init))
def beta_init(m, k, i):
beta = np.array([[-10, -0.75, -2],[0, 0, -1]])
return beta[k, i]
m.del_component('beta')
m.add_component('beta', pe.Param(range(N), range(n), mutable=True, initialize=beta_init))
def v_init(m, i):
v = np.ones((1,n))
return v[0, i]
m.del_component('v')
m.add_component('v', pe.Param(range(n), mutable=True, initialize=v_init)) #
def h_init(m, i):
h = np.ones((1,n))
return h[0, i]
m.del_component('h')
m.add_component('h', pe.Param(range(n), mutable=True, initialize=h_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), mutable=True, initialize=r_init)) # 1 by n
def Epsilon_init(m, i_1, i_2):
Epsilon = np.zeros((n,n));
# last row is all zeros, so stop at n-1
for i in range(n-1):
Epsilon[i, i+1] = 1;
return Epsilon[i_1, i_2]
m.del_component('Epsilon')
m.add_component('Epsilon', pe.Param(range(n), range(n), mutable=True, initialize=Epsilon_init)) # n by n
def gamma_init(m, i):
gamma = np.array([10, 1, 1])
return gamma[i]
m.del_component('gamma')
m.add_component('gamma', pe.Param(range(n), mutable=True, initialize=gamma_init)) # n by n
def w_init(m, i):
w = np.ones((n,1))
return w[i,0]
m.del_component('w')
m.add_component('w', pe.Param(range(n), mutable=True, initialize=w_init)) # n by n
def C_init(m, i):
C = np.ones((n,1))
return C[i, 0]
m.del_component('C')
m.add_component('C', pe.Param(range(n), mutable=True, initialize=C_init)) # n by n
def R_matrix_init(m, k, i):
R_matrix = np.array([[0,1,0], [0,0,1], [0,0,0]])
return R_matrix[k, i]
m.del_component('R_matrix')
m.add_component('R_matrix', pe.Param(range(n), range(n), mutable=True, initialize=R_matrix_init))
return m
def add_inputs(m, Q, q):
n = pe.value(m.n)
N = pe.value(m.N)
m.partition_dim_dict = {'0':6}
m.del_component('Q')
m.add_component('Q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=Q))
Q = pe.value(m.Q)
m.del_component('q')
m.add_component('q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=q))
m.del_component('max_dim_S')
m.add_component('max_dim_S', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=6))
# we add Q, l, d, q
# Parameters we need to change if taking inputs from MPZP: Q, l, d, q
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))
max_dim_S = m.max_dim_S.value
def A_0_init(m, s, i):
A_0 = np.array([[ 1., 0., 0.], [-1., 0., 0.], [ 0., -1., 0.], [ 0., 0., -1.], [ 0., 1., 0.], [ 0., 0., 1.]])
return A_0[s, i]
m.del_component('A_0')
m.add_component('A_0', pe.Param(range(max_dim_S), range(n), mutable=True, initialize=A_0_init))
def b_0_init(m, s):
b_0 = np.array([[ 1.0e+04], [-1.0e+00], [-1.0e+00], [-1.0e+00], [ 9.9e+01], [ 9.9e+01]])
return b_0[s]
m.del_component('b_0')
m.add_component('b_0', pe.Param(range(max_dim_S), mutable=True, initialize=b_0_init))
def l_init(m, ix):
l = np.array([0.0,0.0,0.0,1.0,0.0,0.0])
return l[ix]
m.del_component('l')
m.add_component('l', pe.Param(range(n*Q*N), mutable=True, within=pe.Reals, initialize=l_init))
def d_q_s_init(m, q, s):
return 1
m.del_component('d_q_s')
m.add_component('d_q_s', pe.Param(range(Q), range(max_dim_S), mutable=True, initialize=d_q_s_init))
return m
Q = 1
q = 1
m = setup_model()
m = add_inputs(m, Q, q)
n=pe.value(m.n)
Q=pe.value(m.Q)
N=pe.value(m.N)
m.del_component('Z')
m.add_component('Z', pe.Var(range(n), range(n*Q*N), within=pe.Reals))
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))
import pyomo.environ as pe
from enum import Enum
import numpy as np
import utils
import scipy.optimize as sp_opt
import cloudpickle
import pickle
import pyomo.contrib.pyros as pyros
import generate_constraints_clcctm as gen_constr
import generate_objective_function_clcctm as gen_obj
from pyomo.core.expr import current as EXPR
from pyomo.core.base.constraint import ConstraintList
class Geometry(Enum):
'''
Enum defining uncertainty set geometries
'''
ZINEAR = 1
CONVEX_NONZINEAR = 2
GENERAZ_NONZINEAR = 3
DISCRETE_SCENARIOS = 4
class CustomSet(pyros.UncertaintySet):
def __init__(self, m):
# "m" must be a pyomo model object
self.m = m
# bounds on y
self.type = 'CustomSet'
bounds_d_q_s = zip([0]*m.Q.value*m.max_dim_S.value, [1]*m.Q.value*m.max_dim_S.value)
bounds_d = zip([0]*m.Q.value, [1]*m.Q.value)
m.bounds_y = [(0,10), (0,1), (0,1)]
# === Construct the desirable uncertainty set ===
self.bounds = list(bounds_d_q_s) + list(bounds_d) + m.bounds_y
@property
def dim(self):
m = self.m
return m.Q.value*m.max_dim_S.value + m.Q.value + m.n.value
@property
def geometry(self):
return Geometry.GENERAZ_NONZINEAR
@property
def parameter_bounds(self):
return self.bounds
def set_as_constraint(self, uncertain_params):
m = self.m
## separate the uncertain parameters for convenience
y_start_ix = m.Q.value*m.max_dim_S.value + m.Q.value
d_start_ix = m.Q.value*m.max_dim_S.value
d_q_s_start_ix = 0
conlist = ConstraintList()
conlist.construct()
Q = m.Q.value
max_dim_S = m.max_dim_S.value
n = pe.value(m.n)
N = pe.value(m.N)
epsilon = 0
# linear inequalities
for q in range(Q):
A = eval('m.A_'+str(q))
b = eval('m.b_'+str(q))
S_q = m.partition_dim_dict[str(q)]
for s in range(S_q):
a_s = A[s,:]
a_s_sum = 0
for j in range(n):
new_term = pe.prod([A[s, j], uncertain_params[y_start_ix+j]])
a_s_sum = pe.quicksum([a_s_sum, new_term])
b_s = b[s]
bounds = list(m.bounds_y)
c = list(a_s.value)
m_min = -10000000
d_q_s_ix = q*max_dim_S + s
constraint = a_s_sum - b_s >= epsilon + (m_min - epsilon)*uncertain_params[d_q_s_start_ix+d_q_s_ix]
conlist.add(constraint)
# constraints for extra d_q_s parameters
for q in range(Q):
S_q = m.partition_dim_dict[str(q)]
if S_q < max_dim_S:
# compute difference to find the extra d_q_s parameters
diff = max_dim_S - S_q
for s in range(S_q, max_dim_S):
d_q_s_ix = q*max_dim_S + s
conlist.add(uncertain_params[d_q_s_ix] == 1)
# product constraint (Q of these)
constraint_tol = 0.0001
for q in range(Q):
S_q = m.partition_dim_dict[str(q)]
prod_expression = 1
for s in range(S_q):
d_q_s_ix = q*max_dim_S + s
prod_expression = pe.prod([prod_expression, uncertain_params[d_q_s_ix]])
conlist.add(prod_expression <= uncertain_params[d_start_ix+q] + constraint_tol)
conlist.add(prod_expression >= uncertain_params[d_start_ix+q] - constraint_tol)
# Total number of constraints = max_dim_S*Q + Q, because if a linear inequality was not created, then a constraint should have been created for the extra d_q_s parameter
return conlist
def create_polyhedral_unc_set(m):
n = m.n.value
N = m.N.value
Q = m.Q.value
max_dim_S = m.max_dim_S.value
A_col_dim = max_dim_S*Q + Q + n
A = np.zeros((0, A_col_dim))
b = np.zeros((0, 1))
epsilon = 0.0001
# Add exclusion constraint row first
exclusion_constraint_row_A_lt = np.concatenate([np.zeros((1, max_dim_S*Q)), np.ones((1, Q)), np.zeros((1, n))], axis=1)
exclusion_constraint_row_b_lt = np.array([[1+epsilon]]) # This epsilon is necessary regardless of the epsilon flag to avoid numerical issues with satisfying this constraint
A = np.concatenate([A, exclusion_constraint_row_A_lt], axis=0)
b = np.concatenate([b, exclusion_constraint_row_b_lt], axis=0)
# Add the other inequality
exclusion_constraint_row_A_gt = np.concatenate([np.zeros((1, max_dim_S*Q)), -1*np.ones((1, Q)), np.zeros((1, n))], axis=1)
exclusion_constraint_row_b_gt = np.array([[-1+epsilon]]) # This epsilon is necessary regardless of the epsilon flag to avoid numerical issues with satisfying this constraint
A = np.concatenate([A, exclusion_constraint_row_A_gt], axis=0)
b = np.concatenate([b, exclusion_constraint_row_b_gt], axis=0)
for q in range(Q):
# upper bound
exclusion_constraint_row_A = np.zeros((1, max_dim_S*Q + Q + n))
exclusion_constraint_row_A[0, max_dim_S*Q +q] = 1
exclusion_constraint_row_b = np.array([[1]])
A = np.concatenate([A, exclusion_constraint_row_A], axis=0)
b = np.concatenate([b, exclusion_constraint_row_b], axis=0)
# lower bound
exclusion_constraint_row_A = np.zeros((1, max_dim_S*Q + Q + n))
exclusion_constraint_row_A[0, max_dim_S*Q + q] = -1
exclusion_constraint_row_b = np.array([[0]]) # This is intended to prevent delta_q = 0, because that negates the constraints.
A = np.concatenate([A, exclusion_constraint_row_A], axis=0)
b = np.concatenate([b, exclusion_constraint_row_b], axis=0)
# uncertainty set for y
gamma = [m.gamma[i].value for i in m.gamma._index_set]
h = [m.h[i].value for i in m.h._index_set]
max_densities = np.array([gamma_h[0]*gamma_h[1] for gamma_h in zip(gamma, h)])
nominal_y = 0.5
for i in range(n):
# upper bound
exclusion_constraint_row_A = np.zeros((1, max_dim_S*Q + Q + n))
exclusion_constraint_row_A[0, max_dim_S*Q + Q + i] = 1
exclusion_constraint_row_b = np.array([[max_densities[i]]])
A = np.concatenate([A, exclusion_constraint_row_A], axis=0)
b = np.concatenate([b, exclusion_constraint_row_b], axis=0)
# lower bound
exclusion_constraint_row_A = np.zeros((1, max_dim_S*Q + Q + n))
exclusion_constraint_row_A[0, max_dim_S*Q + Q + i] = -1
exclusion_constraint_row_b = np.array([[0]])
A = np.concatenate([A, exclusion_constraint_row_A], axis=0)
b = np.concatenate([b, exclusion_constraint_row_b], axis=0)
Q = m.Q.value
max_dim_S = m.max_dim_S.value
n = pe.value(m.n)
N = pe.value(m.N)
for q in range(Q):
for s in range(max_dim_S):
# upper bound
new_row_A = np.zeros((1, max_dim_S*Q + Q + n))
new_row_A[0, (max_dim_S*q) + s] = 1
new_b = np.array([[1]])
A = np.concatenate([A, new_row_A], axis=0)
b = np.concatenate([b, new_b], axis=0)
# lower bound
new_row_A = np.zeros((1, max_dim_S*Q + Q + n))
new_row_A[0, (max_dim_S*q) + s] = -1
new_b = np.array([[0]])
A = np.concatenate([A, new_row_A], axis=0)
b = np.concatenate([b, new_b], axis=0)
lhs_coefficient_mat = A
rhs_vec = b.flatten().tolist()
return lhs_coefficient_mat, rhs_vec
# test the intersection alone
lhs_coefficients_mat, rhs_vec = create_polyhedral_unc_set(m)
polyhedral_set = pyros.PolyhedralSet(lhs_coefficients_mat, rhs_vec)
custom_set = CustomSet(m)
uncertainty_set = pyros.IntersectionSet(polyhedral_set = polyhedral_set, custom_set= custom_set)
d_q_s_values = utils.extract_values_into_array(m.d_q_s)
d_values = utils.extract_values_into_array(m.d)
yues = utils.extract_values_into_array(m.y)
unc_param_values = d_q_s_values+d_values+yues
# Check if nominal point is in the intersection
if uncertainty_set.point_in_set(unc_param_values):
print('nominal point is in the intersection!')
uncertain_flag = True
MA_27_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_q_s, m.d, m.y]
lhs_coefficients_mat, rhs_vec = create_polyhedral_unc_set(m)
polyhedral_set = pyros.PolyhedralSet(lhs_coefficients_mat, rhs_vec)
partition_logic_set = CustomSet(m)
uncertainty_set = pyros.IntersectionSet(polyhedral_set = polyhedral_set, partition_logic_set = partition_logic_set)
# === Designate local and global NLP solvers ===
if MA_27_flag:
ipopt_solver = pe.SolverFactory('ipopt')
ipopt_solver.options['OF_hsllib'] = 'libcoinhsl.dylib'
ipopt_solver.options['OF_linear_solver'] = 'ma27'
else:
ipopt_solver = pe.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.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":True,
"bypass_global_separation":False
})
else:
if MA_27_flag:
solver = pe.SolverFactory('ipopt')
solver.options['OF_hsllib'] = 'libcoinhsl.dylib'
solver.options['OF_linear_solver'] = 'ma27'
else:
solver = pe.SolverFactory('ipopt')
results = solver.solve(m, tee=True)
