jump-dev / HiGHS.jl

A Julia interface to the HiGHS solver
https://highs.dev
MIT License
108 stars 15 forks source link

Duals from QP can be incorrect #157

Closed odow closed 1 year ago

odow commented 1 year ago
julia> model = HiGHS.Optimizer()
A HiGHS model with 0 columns and 0 rows.

julia> x, y = MOI.add_variables(model, 2)
Running HiGHS 1.5.1 [date: 1970-01-01, git hash: 93f1876e4]
Copyright (c) 2023 HiGHS under MIT licence terms
2-element Vector{MathOptInterface.VariableIndex}:
 MOI.VariableIndex(1)
 MOI.VariableIndex(2)

julia> MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

julia> f = 2.0*x*x + 1.0*x*y + 1.0*y*y + x + y
0.0 + 1.0 MOI.VariableIndex(1) + 1.0 MOI.VariableIndex(2) + 2.0 MOI.VariableIndex(1)² + 1.0 MOI.VariableIndex(1)*MOI.VariableIndex(2) + 1.0 MOI.VariableIndex(2)²

julia> MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)

julia> c = MOI.add_constraint(model, 1.0 * x + y, MOI.LessThan(-1.0))
MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}(1)

julia> MOI.optimize!(model)
0, 1.000000, 1, 0.000400, 0.000000, 0, 0.000000, 0.000000
3, -0.125000, 1, 0.000482, 0.000000, 0, 0.000000, 1.000000
Model   status      : Optimal
QP ASM    iterations: 3
Objective value     : -1.2500000000e-01
HiGHS run time      :          0.00

julia> MOI.get(model, MOI.ConstraintDual(), c)
0.0

The last dual is incorrect. It's because the basis status is incorrectly reported and we return 0.0 instead to the true dual (HiGHS does find -0.75 for rowdual).

julia> model.solution.rowdual
1-element Vector{Float64}:
 -0.7500000624999998

julia> model.solution.rowstatus
1-element Vector{Int32}:
 0

x-ref https://github.com/ERGO-Code/HiGHS/issues/1249

odow commented 1 year ago

With

diff --git a/src/interfaces/highs_c_api.cpp b/src/interfaces/highs_c_api.cpp
index 821b39e18..5806b60f1 100644
--- a/src/interfaces/highs_c_api.cpp
+++ b/src/interfaces/highs_c_api.cpp
@@ -463,6 +463,9 @@ HighsInt Highs_getBasis(const void* highs, HighsInt* col_status,
   for (HighsInt i = 0; i < (HighsInt)basis.row_status.size(); i++) {
     row_status[i] = (HighsInt)basis.row_status[i];
   }
+  if (!basis.valid) {
+    return kHighsStatusError;
+  }
   return kHighsStatusOk;
 }

and

#include "interfaces/highs_c_api.h"

#include <stdlib.h>
#include <assert.h>

int main(int argc, char *argv[]) {
    void* highs = Highs_create();
    double inf = Highs_getInfinity(highs);
    int num_col = 2;
    int num_row = 1;
    int num_nz = 2;
    int q_num_nz = 3;
    int a_format = kHighsMatrixFormatColwise;
    int q_format = kHighsHessianFormatTriangular;
    int sense = kHighsObjSenseMinimize;
    double offset = 0.0;
    double col_cost[] = {1.0, 1.0};
    double col_lower[] = {-inf, -inf};
    double col_upper[] = {inf, inf};
    double row_lower[] = {-inf};
    double row_upper[] = {-1.0};
    int a_start[] = {0, 1};
    int a_index[] = {0, 0};
    double a_value[] = {1.0, 1.0};
    int q_start[] = {0, 2};
    int q_index[] = {0, 1, 1};
    double q_value[] = {4.0, 1.0, 2.0};
    int integrality[] = {kHighsVarTypeContinuous, kHighsVarTypeContinuous};
    int ret;
    ret = Highs_passModel(highs, num_col, num_row, num_nz, q_num_nz, a_format,
                          q_format, sense, offset, col_cost, col_lower,
                          col_upper, row_lower, row_upper, a_start, a_index,
                          a_value, q_start, q_index, q_value, NULL);
    assert(ret == 0);
    ret = Highs_run(highs);
    assert(ret == 0);
    int col_basis_status[] = {-123, -123};
    int row_basis_status[] = {-123};
    ret = Highs_getBasis(highs, col_basis_status, row_basis_status);
    assert(ret == 0);
    Highs_destroy(highs);
    return 0;
}

I get

(base) oscar@Oscars-MBP examples % export highs=/Users/Oscar/Documents/Code/HiGHS/build;
(base) oscar@Oscars-MBP examples % export DYLD_LIBRARY_PATH=$highs/lib; ./main
Running HiGHS 1.5.1 [date: 2023-04-12, git hash: 46574b917-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
0, 1.000000, 1, 0.000768, 0.000000, 0, 0.000000, 0.000000
3, -0.125000, 1, 0.000844, 0.000000, 0, 0.000000, 1.000000
Model   status      : Optimal
QP ASM    iterations: 3
Objective value     : -1.2500000000e-01
HiGHS run time      :          0.00
Assertion failed: (ret == 0), function main, file main.c, line 40.
zsh: abort      ./main

Let me know if you want me to make a PR. At the very least, we should signal that the basis isn't valid instead of always returning kHighsStatusOk.