PyROS doesn't find feasible solution, when clearly one exists
57 views
Skip to first unread message
huder
Jul 24, 2022, 11:09:10 PM7/24/22
to Pyomo Developers
I’m running a fairly simple robust optimization problem, to which I know there exists a feasible solution. The separation problem seems to find an optimal solution many times, but for some reason when it reports its results back to the master problem, it’s not able to find a feasible solution.
I changed the solver settings to make the problem easier: “objective_focus = pyros.ObjectiveType.nominal” and “solve_master_globally: False”, and objective = 0, but to no avail. I know that the nominal values for the uncertain parameters are in the uncertainty set, or else PyROS would complain.
I reduced the problem dimensions as much as possible to try to figure out why the solver can’t find this solution. In the code excerpt below, I show that there exists a feasible solution by showing that the constraints are satisfied, but “appsi_ipopt” is unable to find this solution.
Is there anything else I can do to help this problem along? Or is this what I should expect when using IPOPT as the local and global solver using PyROS for a problem like this?
Thanks.
import pyomo.environ as pe
import numpy as np
def create_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):
x_0_val = 0.5*np.ones((1,3))
return x_0_val[0,i]
m.del_component('x_0_val')
m.add_component('x_0_val', pe.Param(range(n), mutable=True, initialize=x_0_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)) #
m.del_component('T_s')
m.add_component('T_s', pe.Param(mutable=True, within=pe.Reals, initialize=1)) #
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 gamma_init(m, i):
gamma = np.array([1000, 100, 100])
return gamma[i]
m.del_component('gamma')
m.add_component('gamma', pe.Param(range(n), mutable=True, initialize=gamma_init)) # n by n
return m
m=create_model()
q = 1
Q = 2
n = pe.value(m.n)
N = pe.value(m.N)
m.del_component('q')
m.add_component('q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=q))
m.del_component('Q')
m.add_component('Q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=Q))
def d_init(m, ix):
if ix == m.q.value:
return 1
return 0
m.del_component('d')
m.add_component('d', pe.Param(range(m.Q.value), mutable=True, initialize=d_init))
l_optimal = np.array([[0,0,0,1,0,0]])
l_value = np.concatenate([np.zeros((1, N*n*(q-1))), l_optimal, np.zeros((1, (Q-q)*n*N))], axis=1)
def l_init(m, ix):
return l_value[0, ix]
m.del_component('l')
m.add_component('l', pe.Param(range(n*Q*N), mutable=True, initialize=l_init))
m.del_component('L')
m.add_component('L', pe.Var(range(n), range(n*Q*N), within=pe.Reals))
n = pe.value(m.n)
N = pe.value(m.N)
L_1 = np.eye(n,n)
L_2 = np.zeros((n,n))
L_2[1, 0] = 1
L_2[2, 1] = 1
L_optimal = np.concatenate([L_1, L_2], axis=1)
L_value = np.concatenate([np.zeros((n, q*n*N)), L_optimal, np.zeros((n, (Q-q-1)*n*N))], axis=1)
def assign_values(model, L_value):
[m, n] = np.shape(L_value)
for i in range(m):
for j in range(n):
model.L[i,j].value = L_value[i,j]
return
assign_values(m, L_value)
m.del_component('c1')
m.add_component('c1', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.d[0]*m.L[0,2]*m.x_0_val[2] + m.d[1]*m.L[0,6]*m.x_0_val[0] + m.d[1]*m.L[0,7]*m.x_0_val[1] + m.d[1]*m.L[0,8]*m.x_0_val[2] + (m.d[0]*m.l[0] + m.d[1]*m.l[6])))
m.del_component('c2')
m.add_component('c2', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[0,3] + m.x_0_val[1]*m.d[0]*m.L[0,4] + m.d[0]*m.L[0,5]*m.x_0_val[2] + m.d[1]*m.L[0,9]*m.x_0_val[0] + m.d[1]*m.L[0,10]*m.x_0_val[1] + m.d[1]*m.L[0,11]*m.x_0_val[2] + (m.d[0]*m.l[3] + m.d[1]*m.l[9])))
m.del_component('c3')
m.add_component('c3', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.d[0]*m.L[1,2]*m.x_0_val[2] + m.d[1]*m.L[1,6]*m.x_0_val[0] + m.d[1]*m.L[1,7]*m.x_0_val[1] + m.d[1]*m.L[1,8]*m.x_0_val[2] + (m.d[0]*m.l[1] + m.d[1]*m.l[7])))
m.del_component('c4')
m.add_component('c4', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[1,3] + m.x_0_val[1]*m.d[0]*m.L[1,4] + m.d[0]*m.L[1,5]*m.x_0_val[2] + m.d[1]*m.L[1,9]*m.x_0_val[0] + m.d[1]*m.L[1,10]*m.x_0_val[1] + m.d[1]*m.L[1,11]*m.x_0_val[2] + (m.d[0]*m.l[4] + m.d[1]*m.l[10])))
m.del_component('c5')
m.add_component('c5', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[2,0] + m.x_0_val[1]*m.d[0]*m.L[2,1] + m.d[0]*m.L[2,2]*m.x_0_val[2] + m.d[1]*m.L[2,6]*m.x_0_val[0] + m.d[1]*m.L[2,7]*m.x_0_val[1] + m.d[1]*m.L[2,8]*m.x_0_val[2] + (m.d[0]*m.l[2] + m.d[1]*m.l[8])))
m.del_component('c6')
m.add_component('c6', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[2,3] + m.x_0_val[1]*m.d[0]*m.L[2,4] + m.d[0]*m.L[2,5]*m.x_0_val[2] + m.d[1]*m.L[2,9]*m.x_0_val[0] + m.d[1]*m.L[2,10]*m.x_0_val[1] + m.d[1]*m.L[2,11]*m.x_0_val[2] + (m.d[0]*m.l[5] + m.d[1]*m.l[11])))
m.del_component('c7')
m.add_component('c7', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.d[0]*m.L[0,2]*m.x_0_val[2] + m.d[1]*m.L[0,6]*m.x_0_val[0] + m.d[1]*m.L[0,7]*m.x_0_val[1] + m.d[1]*m.L[0,8]*m.x_0_val[2] + (m.d[0]*m.l[0] + m.d[1]*m.l[6]) <= m.x_0_val[0]))
m.del_component('c8')
m.add_component('c8', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[0,3] + m.x_0_val[1]*m.d[0]*m.L[0,4] + m.d[0]*m.L[0,5]*m.x_0_val[2] + m.d[1]*m.L[0,9]*m.x_0_val[0] + m.d[1]*m.L[0,10]*m.x_0_val[1] + m.d[1]*m.L[0,11]*m.x_0_val[2] + (m.d[0]*m.l[3] + m.d[1]*m.l[9]) <= m.x_0_val[0] + m.r[0]*m.T_s - m.T_s*(m.d[0]*m.l[0] + m.d[1]*m.l[6] + m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.x_0_val[2]*m.d[0]*m.L[0,2] + m.x_0_val[0]*m.d[1]*m.L[0,6] + m.x_0_val[1]*m.d[1]*m.L[0,7] + m.x_0_val[2]*m.d[1]*m.L[0,8])))
m.del_component('c8')
m.add_component('c8', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[0,3] + m.x_0_val[1]*m.d[0]*m.L[0,4] + m.d[0]*m.L[0,5]*m.x_0_val[2] + m.d[1]*m.L[0,9]*m.x_0_val[0] + m.d[1]*m.L[0,10]*m.x_0_val[1] + m.d[1]*m.L[0,11]*m.x_0_val[2] + (m.d[0]*m.l[3] + m.d[1]*m.l[9]) <= m.x_0_val[0] + m.r[0]*m.T_s - m.T_s*(m.d[0]*m.l[0] + m.d[1]*m.l[6] + m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.x_0_val[2]*m.d[0]*m.L[0,2] + m.x_0_val[0]*m.d[1]*m.L[0,6] + m.x_0_val[1]*m.d[1]*m.L[0,7] + m.x_0_val[2]*m.d[1]*m.L[0,8])))
m.del_component('c9')
m.add_component('c9', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.d[0]*m.L[1,2]*m.x_0_val[2] + m.d[1]*m.L[1,6]*m.x_0_val[0] + m.d[1]*m.L[1,7]*m.x_0_val[1] + m.d[1]*m.L[1,8]*m.x_0_val[2] + (m.d[0]*m.l[1] + m.d[1]*m.l[7]) <= m.x_0_val[1]))
m.del_component('c10')
m.add_component('c10', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[1,3] + m.x_0_val[1]*m.d[0]*m.L[1,4] + m.d[0]*m.L[1,5]*m.x_0_val[2] + m.d[1]*m.L[1,9]*m.x_0_val[0] + m.d[1]*m.L[1,10]*m.x_0_val[1] + m.d[1]*m.L[1,11]*m.x_0_val[2] + (m.d[0]*m.l[4] + m.d[1]*m.l[10]) <= m.x_0_val[1] + m.v[1]/m.h[1]*(1 - m.r[1])*(m.v[1]/m.h[1])*m.T_s*(m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.d[0]*m.L[0,2]*m.x_0_val[2] + m.d[1]*m.L[0,6]*m.x_0_val[0] + m.d[1]*m.L[0,7]*m.x_0_val[1] + m.d[1]*m.L[0,8]*m.x_0_val[2] + (m.d[0]*m.l[0] + m.d[1]*m.l[6])) - m.T_s*(m.d[0]*m.l[1] + m.d[1]*m.l[7] + m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.x_0_val[2]*m.d[0]*m.L[1,2] + m.x_0_val[0]*m.d[1]*m.L[1,6] + m.x_0_val[1]*m.d[1]*m.L[1,7] + m.x_0_val[2]*m.d[1]*m.L[1,8])))
m.del_component('c11')
m.add_component('c11', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[2,0] + m.x_0_val[1]*m.d[0]*m.L[2,1] + m.d[0]*m.L[2,2]*m.x_0_val[2] + m.d[1]*m.L[2,6]*m.x_0_val[0] + m.d[1]*m.L[2,7]*m.x_0_val[1] + m.d[1]*m.L[2,8]*m.x_0_val[2] + (m.d[0]*m.l[2] + m.d[1]*m.l[8]) <= m.x_0_val[2]))
m.del_component('c12')
m.add_component('c12', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[2,3] + m.x_0_val[1]*m.d[0]*m.L[2,4] + m.d[0]*m.L[2,5]*m.x_0_val[2] + m.d[1]*m.L[2,9]*m.x_0_val[0] + m.d[1]*m.L[2,10]*m.x_0_val[1] + m.d[1]*m.L[2,11]*m.x_0_val[2] + (m.d[0]*m.l[5] + m.d[1]*m.l[11]) <= m.x_0_val[2] + m.v[2]/m.h[2]*(1 - m.r[2])*(m.v[2]/m.h[2])*m.T_s*(m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.d[0]*m.L[1,2]*m.x_0_val[2] + m.d[1]*m.L[1,6]*m.x_0_val[0] + m.d[1]*m.L[1,7]*m.x_0_val[1] + m.d[1]*m.L[1,8]*m.x_0_val[2] + (m.d[0]*m.l[1] + m.d[1]*m.l[7])) - m.T_s*(m.d[0]*m.l[2] + m.d[1]*m.l[8] + m.x_0_val[0]*m.d[0]*m.L[2,0] + m.x_0_val[1]*m.d[0]*m.L[2,1] + m.x_0_val[2]*m.d[0]*m.L[2,2] + m.x_0_val[0]*m.d[1]*m.L[2,6] + m.x_0_val[1]*m.d[1]*m.L[2,7] + m.x_0_val[2]*m.d[1]*m.L[2,8])))
# ===== show that constraints are satisfied for the above values for the decision variable, and fixed parameter "l"
def all_satisfied():
Truth_array = []
for i in range(n):
for k in range(N):
eval_str = 'Truth_array.append(pe.value(m.c'+str((i*2)+k+1)+'.expr))'
eval(eval_str)
if all(Truth_array):
return True
else:
return False
if all_satisfied():
print('feasible')
else:
print('infeasible')
# === Specify the objective function ===
m.del_component('objective_function')
m.add_component('objective_fn_expression', pe.Objective(expr = 0, sense = pe.minimize))
import pyomo.contrib.pyros as pyros
# === Instantiate the PyROS solver object ===
pyros_solver = pe.SolverFactory("pyros")
gamma = [pe.value(m.gamma[0]), pe.value(m.gamma[1]), pe.value(m.gamma[2])]
h = [pe.value(m.h[0]), pe.value(m.h[1]), pe.value(m.h[2])]
max_densities = np.array([gamma_h[0]*gamma_h[1] for gamma_h in zip(gamma, h)])
bounds_x_0_val = list(zip([0]*m.n.value, max_densities))
bounds_d = [(0,1)]*m.Q.value
# === Specify which parameters are uncertain ===
uncertain_parameters = [m.d, m.x_0_val]
# === Construct the desirable uncertainty set ===
uncertainty_set = pyros.BoxSet(bounds = bounds_x_0_val + bounds_d)
# === 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.L]
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.nominal,
"solve_master_globally": False,
"load_solution":False
})
# === Query results ===
time = results.time
iterations = results.iterations
termination_condition = results.pyros_termination_condition
objective = results.final_objective_value
# === Print some results ===
single_stage_final_objective = round(objective, -1)
print("Final objective value: %s" % single_stage_final_objective)
import numpy as np
def create_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):
x_0_val = 0.5*np.ones((1,3))
return x_0_val[0,i]
m.del_component('x_0_val')
m.add_component('x_0_val', pe.Param(range(n), mutable=True, initialize=x_0_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)) #
m.del_component('T_s')
m.add_component('T_s', pe.Param(mutable=True, within=pe.Reals, initialize=1)) #
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 gamma_init(m, i):
gamma = np.array([1000, 100, 100])
return gamma[i]
m.del_component('gamma')
m.add_component('gamma', pe.Param(range(n), mutable=True, initialize=gamma_init)) # n by n
return m
m=create_model()
q = 1
Q = 2
n = pe.value(m.n)
N = pe.value(m.N)
m.del_component('q')
m.add_component('q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=q))
m.del_component('Q')
m.add_component('Q', pe.Param(mutable=True, within=pe.NonNegativeIntegers, initialize=Q))
def d_init(m, ix):
if ix == m.q.value:
return 1
return 0
m.del_component('d')
m.add_component('d', pe.Param(range(m.Q.value), mutable=True, initialize=d_init))
l_optimal = np.array([[0,0,0,1,0,0]])
l_value = np.concatenate([np.zeros((1, N*n*(q-1))), l_optimal, np.zeros((1, (Q-q)*n*N))], axis=1)
def l_init(m, ix):
return l_value[0, ix]
m.del_component('l')
m.add_component('l', pe.Param(range(n*Q*N), mutable=True, initialize=l_init))
m.del_component('L')
m.add_component('L', pe.Var(range(n), range(n*Q*N), within=pe.Reals))
n = pe.value(m.n)
N = pe.value(m.N)
L_1 = np.eye(n,n)
L_2 = np.zeros((n,n))
L_2[1, 0] = 1
L_2[2, 1] = 1
L_optimal = np.concatenate([L_1, L_2], axis=1)
L_value = np.concatenate([np.zeros((n, q*n*N)), L_optimal, np.zeros((n, (Q-q-1)*n*N))], axis=1)
def assign_values(model, L_value):
[m, n] = np.shape(L_value)
for i in range(m):
for j in range(n):
model.L[i,j].value = L_value[i,j]
return
assign_values(m, L_value)
m.del_component('c1')
m.add_component('c1', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.d[0]*m.L[0,2]*m.x_0_val[2] + m.d[1]*m.L[0,6]*m.x_0_val[0] + m.d[1]*m.L[0,7]*m.x_0_val[1] + m.d[1]*m.L[0,8]*m.x_0_val[2] + (m.d[0]*m.l[0] + m.d[1]*m.l[6])))
m.del_component('c2')
m.add_component('c2', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[0,3] + m.x_0_val[1]*m.d[0]*m.L[0,4] + m.d[0]*m.L[0,5]*m.x_0_val[2] + m.d[1]*m.L[0,9]*m.x_0_val[0] + m.d[1]*m.L[0,10]*m.x_0_val[1] + m.d[1]*m.L[0,11]*m.x_0_val[2] + (m.d[0]*m.l[3] + m.d[1]*m.l[9])))
m.del_component('c3')
m.add_component('c3', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.d[0]*m.L[1,2]*m.x_0_val[2] + m.d[1]*m.L[1,6]*m.x_0_val[0] + m.d[1]*m.L[1,7]*m.x_0_val[1] + m.d[1]*m.L[1,8]*m.x_0_val[2] + (m.d[0]*m.l[1] + m.d[1]*m.l[7])))
m.del_component('c4')
m.add_component('c4', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[1,3] + m.x_0_val[1]*m.d[0]*m.L[1,4] + m.d[0]*m.L[1,5]*m.x_0_val[2] + m.d[1]*m.L[1,9]*m.x_0_val[0] + m.d[1]*m.L[1,10]*m.x_0_val[1] + m.d[1]*m.L[1,11]*m.x_0_val[2] + (m.d[0]*m.l[4] + m.d[1]*m.l[10])))
m.del_component('c5')
m.add_component('c5', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[2,0] + m.x_0_val[1]*m.d[0]*m.L[2,1] + m.d[0]*m.L[2,2]*m.x_0_val[2] + m.d[1]*m.L[2,6]*m.x_0_val[0] + m.d[1]*m.L[2,7]*m.x_0_val[1] + m.d[1]*m.L[2,8]*m.x_0_val[2] + (m.d[0]*m.l[2] + m.d[1]*m.l[8])))
m.del_component('c6')
m.add_component('c6', pe.Constraint(expr = 0 <= m.x_0_val[0]*m.d[0]*m.L[2,3] + m.x_0_val[1]*m.d[0]*m.L[2,4] + m.d[0]*m.L[2,5]*m.x_0_val[2] + m.d[1]*m.L[2,9]*m.x_0_val[0] + m.d[1]*m.L[2,10]*m.x_0_val[1] + m.d[1]*m.L[2,11]*m.x_0_val[2] + (m.d[0]*m.l[5] + m.d[1]*m.l[11])))
m.del_component('c7')
m.add_component('c7', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.d[0]*m.L[0,2]*m.x_0_val[2] + m.d[1]*m.L[0,6]*m.x_0_val[0] + m.d[1]*m.L[0,7]*m.x_0_val[1] + m.d[1]*m.L[0,8]*m.x_0_val[2] + (m.d[0]*m.l[0] + m.d[1]*m.l[6]) <= m.x_0_val[0]))
m.del_component('c8')
m.add_component('c8', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[0,3] + m.x_0_val[1]*m.d[0]*m.L[0,4] + m.d[0]*m.L[0,5]*m.x_0_val[2] + m.d[1]*m.L[0,9]*m.x_0_val[0] + m.d[1]*m.L[0,10]*m.x_0_val[1] + m.d[1]*m.L[0,11]*m.x_0_val[2] + (m.d[0]*m.l[3] + m.d[1]*m.l[9]) <= m.x_0_val[0] + m.r[0]*m.T_s - m.T_s*(m.d[0]*m.l[0] + m.d[1]*m.l[6] + m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.x_0_val[2]*m.d[0]*m.L[0,2] + m.x_0_val[0]*m.d[1]*m.L[0,6] + m.x_0_val[1]*m.d[1]*m.L[0,7] + m.x_0_val[2]*m.d[1]*m.L[0,8])))
m.del_component('c8')
m.add_component('c8', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[0,3] + m.x_0_val[1]*m.d[0]*m.L[0,4] + m.d[0]*m.L[0,5]*m.x_0_val[2] + m.d[1]*m.L[0,9]*m.x_0_val[0] + m.d[1]*m.L[0,10]*m.x_0_val[1] + m.d[1]*m.L[0,11]*m.x_0_val[2] + (m.d[0]*m.l[3] + m.d[1]*m.l[9]) <= m.x_0_val[0] + m.r[0]*m.T_s - m.T_s*(m.d[0]*m.l[0] + m.d[1]*m.l[6] + m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.x_0_val[2]*m.d[0]*m.L[0,2] + m.x_0_val[0]*m.d[1]*m.L[0,6] + m.x_0_val[1]*m.d[1]*m.L[0,7] + m.x_0_val[2]*m.d[1]*m.L[0,8])))
m.del_component('c9')
m.add_component('c9', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.d[0]*m.L[1,2]*m.x_0_val[2] + m.d[1]*m.L[1,6]*m.x_0_val[0] + m.d[1]*m.L[1,7]*m.x_0_val[1] + m.d[1]*m.L[1,8]*m.x_0_val[2] + (m.d[0]*m.l[1] + m.d[1]*m.l[7]) <= m.x_0_val[1]))
m.del_component('c10')
m.add_component('c10', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[1,3] + m.x_0_val[1]*m.d[0]*m.L[1,4] + m.d[0]*m.L[1,5]*m.x_0_val[2] + m.d[1]*m.L[1,9]*m.x_0_val[0] + m.d[1]*m.L[1,10]*m.x_0_val[1] + m.d[1]*m.L[1,11]*m.x_0_val[2] + (m.d[0]*m.l[4] + m.d[1]*m.l[10]) <= m.x_0_val[1] + m.v[1]/m.h[1]*(1 - m.r[1])*(m.v[1]/m.h[1])*m.T_s*(m.x_0_val[0]*m.d[0]*m.L[0,0] + m.x_0_val[1]*m.d[0]*m.L[0,1] + m.d[0]*m.L[0,2]*m.x_0_val[2] + m.d[1]*m.L[0,6]*m.x_0_val[0] + m.d[1]*m.L[0,7]*m.x_0_val[1] + m.d[1]*m.L[0,8]*m.x_0_val[2] + (m.d[0]*m.l[0] + m.d[1]*m.l[6])) - m.T_s*(m.d[0]*m.l[1] + m.d[1]*m.l[7] + m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.x_0_val[2]*m.d[0]*m.L[1,2] + m.x_0_val[0]*m.d[1]*m.L[1,6] + m.x_0_val[1]*m.d[1]*m.L[1,7] + m.x_0_val[2]*m.d[1]*m.L[1,8])))
m.del_component('c11')
m.add_component('c11', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[2,0] + m.x_0_val[1]*m.d[0]*m.L[2,1] + m.d[0]*m.L[2,2]*m.x_0_val[2] + m.d[1]*m.L[2,6]*m.x_0_val[0] + m.d[1]*m.L[2,7]*m.x_0_val[1] + m.d[1]*m.L[2,8]*m.x_0_val[2] + (m.d[0]*m.l[2] + m.d[1]*m.l[8]) <= m.x_0_val[2]))
m.del_component('c12')
m.add_component('c12', pe.Constraint(expr = m.x_0_val[0]*m.d[0]*m.L[2,3] + m.x_0_val[1]*m.d[0]*m.L[2,4] + m.d[0]*m.L[2,5]*m.x_0_val[2] + m.d[1]*m.L[2,9]*m.x_0_val[0] + m.d[1]*m.L[2,10]*m.x_0_val[1] + m.d[1]*m.L[2,11]*m.x_0_val[2] + (m.d[0]*m.l[5] + m.d[1]*m.l[11]) <= m.x_0_val[2] + m.v[2]/m.h[2]*(1 - m.r[2])*(m.v[2]/m.h[2])*m.T_s*(m.x_0_val[0]*m.d[0]*m.L[1,0] + m.x_0_val[1]*m.d[0]*m.L[1,1] + m.d[0]*m.L[1,2]*m.x_0_val[2] + m.d[1]*m.L[1,6]*m.x_0_val[0] + m.d[1]*m.L[1,7]*m.x_0_val[1] + m.d[1]*m.L[1,8]*m.x_0_val[2] + (m.d[0]*m.l[1] + m.d[1]*m.l[7])) - m.T_s*(m.d[0]*m.l[2] + m.d[1]*m.l[8] + m.x_0_val[0]*m.d[0]*m.L[2,0] + m.x_0_val[1]*m.d[0]*m.L[2,1] + m.x_0_val[2]*m.d[0]*m.L[2,2] + m.x_0_val[0]*m.d[1]*m.L[2,6] + m.x_0_val[1]*m.d[1]*m.L[2,7] + m.x_0_val[2]*m.d[1]*m.L[2,8])))
# ===== show that constraints are satisfied for the above values for the decision variable, and fixed parameter "l"
def all_satisfied():
Truth_array = []
for i in range(n):
for k in range(N):
eval_str = 'Truth_array.append(pe.value(m.c'+str((i*2)+k+1)+'.expr))'
eval(eval_str)
if all(Truth_array):
return True
else:
return False
if all_satisfied():
print('feasible')
else:
print('infeasible')
# === Specify the objective function ===
m.del_component('objective_function')
m.add_component('objective_fn_expression', pe.Objective(expr = 0, sense = pe.minimize))
import pyomo.contrib.pyros as pyros
# === Instantiate the PyROS solver object ===
pyros_solver = pe.SolverFactory("pyros")
gamma = [pe.value(m.gamma[0]), pe.value(m.gamma[1]), pe.value(m.gamma[2])]
h = [pe.value(m.h[0]), pe.value(m.h[1]), pe.value(m.h[2])]
max_densities = np.array([gamma_h[0]*gamma_h[1] for gamma_h in zip(gamma, h)])
bounds_x_0_val = list(zip([0]*m.n.value, max_densities))
bounds_d = [(0,1)]*m.Q.value
# === Specify which parameters are uncertain ===
uncertain_parameters = [m.d, m.x_0_val]
# === Construct the desirable uncertainty set ===
uncertainty_set = pyros.BoxSet(bounds = bounds_x_0_val + bounds_d)
# === 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.L]
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.nominal,
"solve_master_globally": False,
"load_solution":False
})
# === Query results ===
time = results.time
iterations = results.iterations
termination_condition = results.pyros_termination_condition
objective = results.final_objective_value
# === Print some results ===
single_stage_final_objective = round(objective, -1)
print("Final objective value: %s" % single_stage_final_objective)
```
Reply all
Reply to author
Forward
0 new messages