7.7 The linear/simplex interface tutorial

In this tutorial we demonstrate the linear/simplex optimization part of API for MATLAB. In most cases linear optimization problems can be specified using the main component of API for MATLAB, the conic optimization model; see Sec. 7.1 (The conic interface tutorial) for an introduction. The separate linear/simplex component documented here is more suitable for users who require:

  • using the simplex algorithm,

  • warm-starting the simplex algorithm,

  • computing a basis solution to a linear problem.

A standard problem in the linear/simplex toolbox is constructed with moseklinmodel and has the form

(7.23)\[\begin{split}\begin{array}{lccccl} \mbox{minimize} & & & c^T x+c^f & & \\ \mbox{subject to} & & & A x & = & b, \\ & b_l & \leq & x & \leq & b_u, \end{array}\end{split}\]

where

  • \(m\) is the number of constraints.

  • \(n\) is the number of decision variables.

  • \(x \in \real^n\) is a vector of decision variables.

  • \(c \in \real^n\) is the linear part of the objective function.

  • \(c^f\in \real\) is a constant term in the objective

  • \(A \in \real^{m \times n}\) is the constraint matrix.

  • \(b \in \real^m\) is the activity of linear constraints.

  • \(b_l \in \real^n\) is the lower bound for the variables.

  • \(b_u \in \real^n\) is the upper bound for the variables.

Lower and upper variable bounds can be infinite, or in other words the corresponding bound may be omitted.

A linear problem in this format can be constructed and input either all at once or incrementally, by adding successive blocks of constraints.

It is also possible to extend the model with new variables, modify bounds, input solutions and status keys. etc. This is a topic of a separate tutorial. TODO

7.7.1 Example LO1

We set up a linear problem in one go. Consider the problem

(7.24)\[\begin{split}\begin{array} {lccccccccl} \mbox{maximize} & 3 x_1 & + & 1 x_2 & + & 5 x_3 & + & 1 x_4 & & \\ \mbox{subject to} & 3 x_1 & + & 1 x_2 & + & 2 x_3 & & & = & 30, \\ & 2 x_1 & + & 1 x_2 & + & 3 x_3 & + & 1 x_4 & \geq & 15, \\ & & & 2 x_2 & & & + & 3 x_4 & \leq & 25, \end{array}\end{split}\]

with variable bounds

\[\begin{split}\begin{array}{ccccc} 0 & \leq & x_1 & \leq & \infty , \\ 0 & \leq & x_2 & \leq & 10, \\ 0 & \leq & x_3 & \leq & \infty ,\\ 0 & \leq & x_4 & \leq & \infty . \end{array}\end{split}\]

We need to reformulate the problem so that the constraints are of the form \(Ax=b\) by introducing nonnegative slack variables:

(7.25)\[\begin{split}\begin{array} {lccccccccccccl} \mbox{maximize} & 3 x_1 & + & 1 x_2 & + & 5 x_3 & + & 1 x_4 & & & & & & \\ \mbox{subject to} & 3 x_1 & + & 1 x_2 & + & 2 x_3 & & & & & & & = & 30, \\ & 2 x_1 & + & 1 x_2 & + & 3 x_3 & + & 1 x_4 & - & s_1 & & & = & 15, \\ & & & 2 x_2 & & & + & 3 x_4 & & & + & s_2 & = & 25, \end{array}\end{split}\]

with variable bounds

\[\begin{split}\begin{array}{ccccc} 0 & \leq & x_1,x_3,x_4,s_1,s_2 & \leq & \infty , \\ 0 & \leq & x_2 & \leq & 10. \end{array}\end{split}\]

We create a model by providing all data at once:

  • the number of variables is 5 and is set as numvar,

  • the objective vector is \(c=[3,1,5,1,0,0]\) and is passed as objective,

  • the objective sense is maximization, and is passed as objsense,

  • the problem name name (optional) is lo1,

  • we also add variable and constraint names for illustration.

Listing 7.15 Setting up a linear model in a single call. Click here to download.
    model = moseklinmodel(...
                  name = "lo1", ...
                  numvar = 6, ...
                  objsense = "maximize", ...
                  objective = [ 3 1 5 1 0 0 ]', ...
                  A = [ 3 1 2 0  0  0 ; ...
                        2 1 3 1 -1  0 ; ...
                        0 2 0 3  0 -1 ], ...
                  b = [ 30 15 25 ]', ...
                  blx = [ 0.0  0.0 0.0 0.0   0.0 -inf ]', ...
                  bux = [ inf 10.0 inf inf   inf  0.0 ]', ...
                  varnames = ["x1" "x2" "x3" "x4" "y1" "y2"]',...
                  conname = ["c1" "c2" "c3"]');

7.7.2 Example with multiple calls (LO2)

We now show how to set up the same model (7.25) adding various elements of the problem in separate calls. For illustration we also split the constraints intwo two blocks and add them separately. The only rule to be observed is that each element can only refer to elements already defined, for example we can not add the objective before defining the number of variables, i.e. the model must remain consistent at all times throughout the process.

Listing 7.16 Setting up a linear model with multiple calls. Click here to download.
    % Create an empty model
    model = moseklinmodel(name = "lo2");

    % Add some variables with bounds
    model.appendvars(6,...
                     bl = [ 0.0  0.0 0.0 0.0   0.0 -inf ]', ...
                     bu = [ inf 10.0 inf inf   inf  0.0 ]');
    model.varname([1:6],["x1" "x2" "x3" "x4" "y1" "y2"]);
    
    % Set the objective
    model.objective("maximize", [ 3 1 5 1 0 0]', objfixterm = 0.0);

    % Constraint 3x1 + x2 + 2x3 == 30
    model.appendcons([ 3 1 2 0 0 0], [30]);

    % The two remaining constraints
    model.appendcons([ 2 1 3 1 -1  0 ; ...
                       0 2 0 3  0 -1 ], ...
                     [ 15 25 ]');

7.7.3 Solving and retrieving the solution

We wrap up with a short demonstration of what to do after the model has been defined, that is how to solve the model and retrieve solutions.

Listing 7.17 Solving and retrieving the solution. Click here to download.
    % Choose the simplex optimizer to solve with
    model.solve(param = [ "MSK_IPAR_OPTIMIZER", "MSK_OPTIMIZER_FREE_SIMPLEX" ]);

    % Check the solution status and fetch the solution
    [hassol, prosta, solsta] = model.hassolution();
    if hassol && solsta == "OPTIMAL"
        fprintf("Solution, problem status = %s, solution status = %s:\n", prosta, solsta);

        % Primal variables
        result = model.getsolution("x");
        x = result(1:4);
        s = result(5:6);
        fprintf("  x[%d] = %.4f\n",[1:4 ; x']);
        fprintf("Slacks:\n");
        fprintf("  s[%d] = %.4f\n",[1:2 ; s']);

        % Other elements of the solution
        fprintf("Duals:\n");
        fprintf("  y[%d] = %.4f\n",[1:3 ; model.getsolution("y")']);
        fprintf("slx[%d] = %.4f\n",[1:6 ; model.getsolution("slx")']);
        fprintf("sux[%d] = %.4f\n",[1:6 ; model.getsolution("sux")']);

        fprintf("Status keys of constraints:\n");
        fprintf("skc[%s] = %s\n",[1:3 ; model.getsolution("skc")']);

        fprintf("Status keys of variables:\n");
        fprintf("skx[%s] = %s\n",[1:6 ; model.getsolution("skx")']);
    end

We note that