scipopt / scip

SCIP - Solving Constraint Integer Programs
Other
392 stars 63 forks source link

Right way of using `SCIPgetDualsolLinear` #49

Closed brmanuel closed 1 year ago

brmanuel commented 1 year ago

I am trying to understand

  1. the right way of using SCIPgetDualsolLinear
  2. the difference between SCIPgetDualsolLinear and SCIPgetDualSolVal

Looking at their documentation, both should do more or less the same

I originally opened a similar issue at PySCIPOpt (issue) they sent me here so I tried to recreate the issue using the C-API:

What did I do

When trying to get the optimal dual values of an LP solution, only SCIPgetDualSolVal seems to return the correct value:

#include <stdio.h>
#include <stdbool.h>

#include "scip/scip.h"
#include <scip/scipdefplugins.h>

int main() {
    // Setup
    SCIP* scip;
    SCIP_CALL_ABORT(SCIPcreate(&scip));
    SCIP_CALL_ABORT(SCIPincludeDefaultPlugins(scip));

    // Turn off features that seem problematic in conjunction with duals:
    // https://www.solveforum.com/forums/threads/solved-how-to-completeley-disable-presolve-in-scip.2697850/
    // https://github.com/scipopt/PySCIPOpt/issues/136
    SCIP_CALL_ABORT(SCIPsetPresolving(scip, SCIP_PARAMSETTING_OFF, true));
    SCIP_CALL_ABORT(SCIPsetHeuristics(scip, SCIP_PARAMSETTING_OFF, true));
    SCIP_CALL_ABORT(SCIPsetIntParam(scip, "propagating/maxrounds", 0));

    SCIP_CALL_ABORT(SCIPcreateProb(scip, "primal", NULL, NULL,
                                 NULL, NULL, NULL, NULL, NULL));
    SCIP_CALL_ABORT(SCIPsetObjsense(scip, SCIP_OBJSENSE_MAXIMIZE));

    // Variables: x,y positive reals
    SCIP_VAR* x;
    SCIP_CALL_ABORT(SCIPcreateVar(scip, &x, "x", 0.0, SCIPinfinity(scip), 1.0,
                                  SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
                                  NULL, NULL, NULL, NULL, NULL));
    SCIP_CALL_ABORT(SCIPaddVar(scip, x));

    SCIP_VAR* y;
    SCIP_CALL_ABORT(SCIPcreateVar(scip, &y, "y", 0.0, SCIPinfinity(scip), 1.0,
                                  SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
                                  NULL, NULL, NULL, NULL, NULL));
    SCIP_CALL_ABORT(SCIPaddVar(scip, y));

    // Constraints
    SCIP_CONS* c1;
    SCIP_VAR* vars[2] = {x, y};
    double coeffs1 [2] = {1.0, 2.0};
    SCIP_CALL_ABORT(SCIPcreateConsLinear(scip, &c1, "c1",
                                         2, vars, coeffs1, 21.0, 21.0, TRUE,
                                         TRUE, TRUE, TRUE, TRUE, FALSE,
                                         FALSE, FALSE, FALSE, FALSE));
    SCIP_CALL_ABORT(SCIPaddCons(scip, c1));

    SCIP_CONS* c2;
    double coeffs2 [2] = {2.0, 1.0};
    SCIP_CALL_ABORT(SCIPcreateConsLinear(scip, &c2, "c2",
                                         2, vars, coeffs2, 24.0, 24.0, TRUE,
                                         TRUE, TRUE, TRUE, TRUE, FALSE,
                                         FALSE, FALSE, FALSE, FALSE));
    SCIP_CALL_ABORT(SCIPaddCons(scip, c2));

    // Solution
    SCIP_CALL_ABORT(SCIPsolve(scip));
    SCIP_SOL* sol = SCIPgetBestSol(scip);
    if (sol == NULL){
        printf("No solution found\n");
    }
    else {

        double x_val = SCIPgetSolVal(scip, sol, x);
        double y_val = SCIPgetSolVal(scip, sol, y);
        double c1_dual = SCIPgetDualsolLinear(scip, c1);
        double c2_dual = SCIPgetDualsolLinear(scip, c2);
        double c1_dual2;
        double c2_dual2;
        SCIP_CALL_ABORT(SCIPgetDualSolVal(scip, c1, &c1_dual2, NULL));
        SCIP_CALL_ABORT(SCIPgetDualSolVal(scip, c2, &c2_dual2, NULL));

        printf("Solution: x = %f, y = %f\n", x_val, y_val);
        printf("getDualsolLinear: c1_dual = %f, c2_dual = %f\n", c1_dual, c2_dual);
        printf("getDualSolVal: c1_dual = %f, c2_dual = %f\n", c1_dual2, c2_dual2);
    }

    // Cleanup
    SCIP_CALL_ABORT(SCIPreleaseVar(scip, &x));
    SCIP_CALL_ABORT(SCIPreleaseVar(scip, &y));
    SCIP_CALL_ABORT(SCIPreleaseCons(scip, &c1));
    SCIP_CALL_ABORT(SCIPreleaseCons(scip, &c2));
    SCIP_CALL_ABORT(SCIPfree(&scip));

    return 0;
}

output:

Solution: x = 9.000000, y = 6.000000
getDualsolLinear: c1_dual = 0.000000, c2_dual = 0.000000
getDualSolVal: c1_dual = 0.333333, c2_dual = 0.333333

The optimal dual values of constraints c1 and c2 are both 1/3, so only getDualSolVal returns the right value but getDualsolLinear doesn't.

Questions

My setup

I'm using SCIP 8.0.3 installed using conda on a Linux machine. More details:

$ scip --version
SCIP version 8.0.3 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 6.0.3] [GitHash: 62fab8a2e3]
Copyright (C) 2002-2022 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)

External libraries: 
  Soplex 6.0.3         Linear Programming Solver developed at Zuse Institute Berlin (soplex.zib.de) [GitHash: f900e3d0]
  CppAD 20180000.0     Algorithmic Differentiation of C++ algorithms developed by B. Bell (github.com/coin-or/CppAD)
  ZLIB 1.2.13          General purpose compression library by J. Gailly and M. Adler (zlib.net)
  GMP 6.2.1            GNU Multiple Precision Arithmetic Library developed by T. Granlund (gmplib.org)
  ZIMPL 3.5.3          Zuse Institute Mathematical Programming Language developed by T. Koch (zimpl.zib.de)
  AMPL/MP 4e2d45c4     AMPL .nl file reader library (github.com/ampl/mp)
  PaPILO 2.1.2         parallel presolve for integer and linear optimization (github.com/scipopt/papilo) [GitHash: 2fe2543]
  bliss 0.77           Computing Graph Automorphism Groups by T. Junttila and P. Kaski (www.tcs.hut.fi/Software/bliss/)
  Ipopt 3.14.12        Interior Point Optimizer developed by A. Waechter et.al. (github.com/coin-or/Ipopt)

Compiler: gcc 12.2.0

Build options:
 ARCH=x86_64
 OSTYPE=Linux-5.15.0-1037-azure
 COMP=GNU 12.2.0
 BUILD=Release
 DEBUGSOL=OFF
 EXPRINT=cppad
 SYM=bliss
 GMP=ON
 IPOPT=ON
 WORHP=OFF
 LPS=spx
 LPSCHECK=OFF
 NOBLKBUFMEM=OFF
 NOBLKMEM=OFF
 NOBUFMEM=OFF
 THREADSAFE=ON;FORCE
 READLINE=OFF
 SANITIZE_ADDRESS=OFF
 SANITIZE_MEMORY=OFF
 SANITIZE_UNDEFINED=OFF
 SANITIZE_THREAD=OFF
 SHARED=ON
 VERSION=8.0.3.0
 API_VERSION=104
 ZIMPL=ON
 ZLIB=ON
mmghannam commented 1 year ago

Hello @brmanuel!

It seems that you have been indirectly introduced to the SCIP's solving stages :)

Before starting the solving process, SCIP transforms the problem into an easier form to work with. That involves creating transformed variables and constraints.

So, here the reason why SCIPgetDualsolLinear returns zeros, is that you pass the original constraints to it which are not part of the transformed problem so they are inactive hence have zero values.

The correct one to use here is SCIPgetDualSolVal as it does in fact call SCIPgetDualsolLinear but calls it on the corresponding transformed constraint.

brmanuel commented 1 year ago

Hi @mmghannam Thank you very much for your answer here and in the related issue! This clears up my confusion.