Adding Indicator constraints for PWL constraints

I have an optimization model of distributed energy system. I want to add a non linear constraint for PEM Elektroliser power calculation depending on hydrogen production. I use pyomo GurobiPersistent Interface for adding PWL constraints. The issue is that depending on a maximum hydrogen production I want to activate certain PWL constraint for corresponding PEM configuration. I am trying to use Indicator constraint for triggering a certain PWL constraint but having continious errors. Is it in general possible to activate certain PWL constraint depending on maximum of another variable?

Here is a snippet of how I tried adding the constraints

# create model
    prob = create_model(data, dt, timesteps, objective, dual=False)
    # prob.write('model.lp', io_options={'symbolic_solver_labels':True})
    prob_filename = os.path.join(result_dir, 'model.lp')
    prob.write(prob_filename) #io_options={'symbolic_solver_labels':True})
    # refresh time stamp string and create filename for logfile
    log_filename = os.path.join(result_dir, '{}.log').format(sce)
    # solve model and read results
    if Solver == "gurobi_persistent":
        # swith to guropi persistent solver for pyomo
        opt = SolverFactory(Solver)
        opt = setup_solver(opt, logfile=log_filename)
        opt.set_instance(prob, symbolic_solver_labels=True)
        opt._solver_model.addGenConstrMax(opt._pyomo_var_to_solver_var_map[prob.max_H2_prod],\
                      [opt._pyomo_var_to_solver_var_map[prob.e_pro_out[tm, 2020, 'CA', 'PEM', 'H2']] for tm in prob.tm])
                   
        # constraint of stack demand from H2 production using Piecewise Linear Approximation
        xpts_stack_1 = []
        ypts_stack_1 = []
        intv_stack_1 = 20.
        xmax_stack_1 = 2575.
        t_stack_1 = 257.
        while t_stack_1 < xmax_stack_1 + intv_stack_1:
            xpts_stack_1.append(t_stack_1)
            ypts_stack_1.append(stack_power(t_stack_1))
            t_stack_1 +=intv_stack_1
     
        for com in prob.com_demand:
            for pro in prob.pro:
                if (com == 'StackDemand' and pro == 'Stack'):
                    for tm in prob.tm:                        
                                opt._solver_model.addGenConstrPWL(opt._pyomo_var_to_solver_var_map[prob.e_pro_out[tm, 2020, 'CA', 'PEM', 'H2']],\
                                    opt._pyomo_var_to_solver_var_map[prob.stack_demand[tm, 2020, 'CA', 'StackDemand', 'Demand']], \
                                        xpts_stack_1, ypts_stack_1, 'pwl_1')                                  
                else:
                    pyomo.Constraint.Skip
        xpts_stack_2 = []
        ypts_stack_2 = []
        intv_stack_2 = 20.
        xmax_stack_2 = 3000.
        t_stack_2 = 300.
        while t_stack_2 < xmax_stack_2 + intv_stack_2:
            xpts_stack_2.append(t_stack_2)
            ypts_stack_2.append(stack_power(t_stack_2))
            t_stack_2 +=intv_stack_2
     
        for com in prob.com_demand:
            for pro in prob.pro:
                if (com == 'StackDemand' and pro == 'Stack'):
                    for tm in prob.tm:                        
                                opt._solver_model.addGenConstrPWL(opt._pyomo_var_to_solver_var_map[prob.e_pro_out[tm, 2020, 'CA', 'PEM', 'H2']],\
                                    opt._pyomo_var_to_solver_var_map[prob.stack_demand[tm, 2020, 'CA', 'StackDemand', 'Demand']], \
                                        xpts_stack_2, ypts_stack_2, 'pwl_2')                                  
                else:
                    pyomo.Constraint.Skip
                                   
        opt._solver_model.addConstr((opt._pyomo_var_to_solver_var_map[prob.max_H2_prod] <= 2600) >> (opt._solver_model.getConstrByName('pwl_1')))
        opt._solver_model.addConstr((opt._pyomo_var_to_solver_var_map[prob.max_H2_prod] >= 2600) >> (opt._solver_model.getConstrByName('pwl_2')))

Jaromił Najman

Gurobi Staff

2023年04月25日 11:34 編集済み

It is not directly possible to model the indicator constraint \(b = 1 \rightarrow y = f(x)\), with \(f\) being a general function, e.g., PWL. However, you can model it with an auxiliary variable \(z\).

You have given \(y = f(x)\) via addGenConstrXxx. You can now model the constraint

\[\begin{align*}
b = 1 \rightarrow z = y\\
\end{align*}\]

and use \(z\) in the rest of your model instead of \(y\). This way, when \(b=1\) then \(z=y\) and it is as if you would have used \(y\), but when \(b=0\) then \(y\) is a free variable because the equality constraint \(z = y\) is not enforced. Thus, the constraint \(y = f(x)\) is redundant because \(y\) is not used anywhere else in the model.

Unfortunately, I am not a Pyomo user, thus I cannot really help on the implementation side but hopefully with the above, you might be able to implement the formulation.

Best regards,
Jaromił

Adding Indicator constraints for PWL constraints

コメント

お探しのものを見つけられませんでしたか？