lp-dual

Models a simple LP and generates its dual. Solves the dual problem with HiGHS and prints, for each primal constraint, its associated dual variable with its value.

The primal LP is taken from "Model building in Mathematical Programming" by H. Paul Williams and is given by

\[ \begin{align*} \min_{x,y} \quad & 3x + 2y \\ \text{s.t.} \quad & x + y \le 4, \\ & 2x + y \le 5, \\ & -x + 4y \ge 2, \\ & x,y \ge 0. \end{align*} \]

The dual problem reads

\[ \begin{align*} \max_{\alpha,\beta,\gamma,\delta,\epsilon} \quad & 4\alpha + 5\beta + 2\gamma \\ \text{s.t.} \quad & \alpha + 2\beta - \gamma + \delta = 3, \\ & \alpha + \beta + 4\gamma + \epsilon = 2, \\ & \alpha, \beta \le 0, \\ & \gamma, \delta, \epsilon \ge 0. \end{align*} \]

#include <iostream>
#include "idol/modeling.h"
#include "idol/mixed-integer/modeling/models/KKT.h"
#include "idol/mixed-integer/optimizers/wrappers/HiGHS/HiGHS.h"
 
using namespace idol;
 
int main(int t_argc, const char** t_argv) {
 
    Env env;
 
    /***************************/
    /* Create the primal model */
    /***************************/
 
    Model primal(env);
 
    // Add primal variables
    const auto x = primal.add_var(0, Inf, Continuous, 3, "x");
    const auto y = primal.add_var(0, Inf, Continuous, 2, "y");
 
    // Add primal constraints
    const auto c1 = primal.add_ctr(x + y <= 4);
    const auto c2 = primal.add_ctr(2 * x + y <= 5);
    const auto c3 = primal.add_ctr(-x + 4 * y >= 2);
 
    /*************************/
    /* Create the dual model */
    /*************************/
 
    Model dual(env);
 
    // Create a KKT reformulator
    // This object can generate various reformulations of the problem such as its dual problem or its KKT conditions.
    Reformulators::KKT reformulator(primal);
 
    // Add the dual problem to the model
    reformulator.add_dual(dual);
 
    // Print both the primal and the dual models
    std::cout << "Primal model:\n" << primal << std::endl;
    std::cout << "Dual model:\n" << dual << std::endl;
 
    /**************************/
    /* Solve the dual problem */
    /**************************/
 
    dual.use(HiGHS());
    dual.optimize();
 
    std::cout << "Dual solution status: " << dual.get_status() << std::endl;
    std::cout << "Dual objective value: " << dual.get_best_obj() << std::endl;
 
    // Get the dual value of each primal constraint
    for (const auto& ctr : primal.ctrs()) {
        const auto& dual_var = reformulator.get_dual_var(ctr);
        std::cout << ctr.name() << " with dual variable " << dual_var.name() << " = " << dual.get_var_primal(dual_var) << std::endl;
    }
 
    // Get the dual value of each primal variable's lower bound constraint
    for (const auto& var : primal.vars()) {
        const auto& dual_var = reformulator.get_dual_var_lb(var);
        std::cout << var.name() << " with dual variable for its lower bound " << dual_var.name() << " = " << dual.get_var_primal(dual_var) << std::endl;
    }
 
    // Primal variables having no upper bound, they do not have a dual variable for their upper bound
 
    return 0;
}