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.

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)