How to help PyROS verify nonempty Intersection Set, when nominal point is clearly in the intersection?

[makansij]

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 nlp in line 1094 of uncertainty_sets.py in 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 

...

ValueError: Solver terminated with an error while checking set intersection non-emptiness.

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 y_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=y_init))

    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 

    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
    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))

from pyomo.core.base.constraint import ConstraintList
import pyomo.environ as pe
from enum import Enum

class Geometry(Enum):
    '''
    Enum defining uncertainty set geometries
    '''
    LINEAR = 1
    CONVEX_NONLINEAR = 2
    GENERAL_NONLINEAR = 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.GENERAL_NONLINEAR

    @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]])  

    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]]) 

    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]])  

        A = np.concatenate([A, exclusion_constraint_row_A], axis=0)
        b = np.concatenate([b, exclusion_constraint_row_b], axis=0)

    # uncertainty set for y
    max_densities = np.array([m.gamma[i].value for i in m.gamma._index_set])

    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 = [m.d_q_s[q_s[0],q_s[1]].value for q_s in m.d_q_s._index_set]
d_values = [m.d[q].value for q in m.d._index_set]
y_values = [m.y[i].value for i in m.y._index_set]
unc_param_values = d_q_s_values+d_values+y_values

# 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 = 0, 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)
    custom_set = CustomSet(m)
    uncertainty_set = pyros.IntersectionSet(polyhedral_set = polyhedral_set, custom_set= custom_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)

[shermanjasonaf]

@makansij I encounter a similar error message:

bad line 121 of /tmp/tmpk5vuidmv.pyomo.nl: 1 array([-10000.])
terminate called after throwing an instance of 'Ipopt::TNLP::INVALID_TNLP'
Ipopt 3.13.2: option_file_name=/tmp/tmpzywh121k_ipopt.opt
ERROR: Solver (ipopt) returned non-zero return code (-6)

This seems to be occuring because the model you use to initialize the CustomSet class has components which are initialized to 1-D Numpy arrays instead of numeric type objects---and these components are carried through to the feasibility problem solved by PyROS to determine whether your IntersectionSet instance represents an empty set. In particular, in your script, you are initializing a Param component 'b_0' with the function b_0_init (in add_inputs), defined through:

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]

This function returns a 1-D Numpy array instead of a numerical value as b_0 is a 2-D Numpy array. Hence, each member of the Param component will be initialized to a 1-D array instead of a float-like value. Was there any particular reason you defined b_0 to be a 2-D array? Alternatively, you could return b_0[s, 0]. If I do that, PyROS now appears to have successfully solved the problem and at a latter stage, raises the exception:

ValueError: No performance constraints identified for the postulated robust optimization problem.

Provided you have resolved your IntersectionSet instance, you may want to check your model constraints.