`linopy` does not pass quadratic objective terms to HiGHS

[fneum]

When I use the quadratic_marginal_cost attribute in PyPSA in combination with HiGHS, I get different results than with Gurobi. HiGHS logs also state that it's solving an LP rather than a QP.

This is because HiGHS does not seem to get the quadratic terms in the objective function (although there seems to be code that should):

https://github.com/PyPSA/linopy/blob/74ce2d2eda87b9a4fe042a6dcc208e2af27b5469/linopy/io.py#L416-L422

[FabianHofmann]

This is weird, as all the tests seem to be fine. Things like

m = linopy.Model()

lower = pd.Series(0, range(10))
x = m.add_variables(lower, name="x")
y = m.add_variables(lower, name="y")

m.add_constraints(x + y >= 10)

m.add_objective(- 2 * x + y + x * x)

m.solve(solver_name="highs")

yield the same result for highs and gurobi. and also the model definition seems totally fine when running test for pypsa networks (however it does not solve with highs). Did it explicitly say that it is trying to solve a LP in your case? and are the dual non-zero?

[fneum]

I think I understand this now.

Let's take a simple LP file with quadratic terms:

min

obj:

-4 x0 -6 x1 + [ + 4 x0 * x0 + 12 x1 * x1 ] / 2

s.t.

c0: +1 x0 +2 x1 >= +4

bounds

+0 <= x0 <= +4
+0 <= x1 <= +4
end

If we run

import highspy
h = highspy.Highs()
h.readModel("test-highs.lp")
h.run()

we get the correct solution:

Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
<HighsStatus.kOk: 0>
0, 16.000001, 0, 0.000281, 0.000000, 0, 0.000000, 0.000000
4, -0.071428, 1, 0.000331, 0.000000, 0, 0.000000, 1.000000
Model   status      : Optimal
QP ASM    iterations: 4
Objective value     : -7.1428571429e-02
HiGHS run time      :          0.00

However, if we specify "ipm" as solver:

import highspy
h = highspy.Highs()
h.readModel("test-highs.lp")
h.setOptionValue("solver", "ipm")
h.run()

we get a different result:

Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
0 rows, 0 cols, 0 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve : Reductions: rows 0(-1); columns 0(-2); elements 0(-2) - Reduced to empty
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Objective value     : -4.0000000000e+01
HiGHS run time      :          0.00

I always assumed that HiGHS would be using the same IPM solver for linear as for quadratic problems, but this does in fact not seem to be the case (https://github.com/ERGO-Code/HiGHS/pull/766). The documentation of HiGHS also notes that implicitly:

solver

Solver option: "simplex", "choose" or "ipm". If "simplex"/"ipm" is chosen then, for a MIP (QP) the integrality constraint (quadratic term) will be ignored
Type: string
Default: "choose"

So, no problem with linopy.

[FabianHofmann]

Great that you cleared that up!