NonConvex + Investment generates unexpected quadratic constraints

Q&A
, ,
1

Hello!

I am currently modeling a small PtX system in oemof.solph with a MeOH synthesis plant as a converter with an investment and nonconvex for the dynamics of the synthesis. When solving with Gurobi, I noticed that the model is reported as having a large number of quadratic constraints:

Optimize a model with 245262 rows, 227748 columns and 604148 nonzeros (Min)
Model fingerprint: 0x041ab51e
Model has 17524 linear objective coefficients
Model has 17516 quadratic constraints

I think the reason for the quadratic constraints is the gradient limits.

Is this behaviour expected for the nonconvex + Investment when using gradient limits?

Cheers,
Anton

meoh_synthesis = solph.components.Converter(
    label="MeOH_Synthesis",
    inputs={
        h2_bus: solph.flows.Flow(),
        el_bus: solph.flows.Flow(),
        co2_bus: solph.flows.Flow(),
        heat_bus: solph.flows.Flow(),
    },
    outputs={
        methanol_bus: solph.flows.Flow(
            min=0.2,  
            max=1.0,   
            nominal_capacity=solph.Investment(
                ep_costs=ep_meoh,  
                maximum=100.0,    
            ),
            nonconvex=solph.NonConvex(
                initial_status=0,
                minimum_uptime=2,
                minimum_downtime=1,
                negative_gradient_limit=0.5,
                positive_gradient_limit=0.5,
            ),

        )
    },
    conversion_factors={
        h2_bus: 1.0 / 0.8,       
        el_bus: 0.04 / 0.8,      
        co2_bus: 211.2 / 0.8,     
        methanol_bus: 1.0,      
        heat_bus: 0.25 / 0.8,     
    },
)
2

Hi @anton !

A straightforward implementation of the combination of Flows with a binary status (historically called NonConvex in solph, although this is technically not always the case) would introduce quadratic (strictly speaking: bilinear) constraints. However, we replace that formulation with a MILP one:

Maybe, Gurobi replaces that because it knows its quadratic solver is faster than working with the binary variables we need. If that’s the case, switching to a different solver (e.g. CBC) should work. If we included quadratic constraints by accident, it would not. Can you check that?

Regards,

Patrik

PS: I know there is an example where a quadratic constraint slipped through a code review. I cannot recall the details, though. Maybe, it is fixed, maybe it is what you describe.