7.1 Linear Optimization

The simplest optimization problem is a purely linear problem. A linear optimization problem (see also Sec. 12.1 (Linear Optimization)) is a problem of the following form:

Minimize or maximize the objective function

$$\sum _{j=0}^{n-1}{c}_j{x}_j+c^f$$

subject to the linear constraints

$$l_k^c\le \sum _{j=0}^{n-1}{a}_{kj}{x}_j\le u_k^c,k=0,\dots ,m-1,$$

and the bounds

$$l_j^x\le x_j\le u_j^x,j=0,\dots ,n-1.$$

The problem description consists of the following elements:

• $m$ and $n$ — the number of constraints and variables, respectively,

• $x$ — the variable vector of length $n$,

• $c$ — the coefficient vector of length $n$

  $$\begin{array}{r}c=\left[\begin{array}{c}{c}_0\\ ⋮\\ c_{n-1}\end{array}\right],\end{array}$$

• $c^f$ — fixed term in the objective,

• $A$ — an $m\times n$ matrix of coefficients

  $$\begin{array}{r}A=\left[\begin{array}{ccc}{a}_{0,0}& \cdots & a_{0,(n-1)}\\ ⋮& \cdots & ⋮\\ a_{(m-1),0}& \cdots & a_{(m-1),(n-1)}\end{array}\right],\end{array}$$

• $l^c$ and $u^c$ — the lower and upper bounds on constraints,

• $l^x$ and $u^x$ — the lower and upper bounds on variables.

Please note that we are using $0$ as the first index: $x_0$ is the first element in variable vector $x$.

The Fusion user does not need to specify all of the above elements explicitly — they will be assembled from the Fusion model.

7.1.1 Example LO1

The following is an example of a small linear optimization problem:

(7.1)$$\begin{array}{r}\begin{array}{lccccccccl}\text{maximize}& 3x_0& +& 1x_1& +& 5x_2& +& 1x_3& & \\ \text{subject to}& 3x_0& +& 1x_1& +& 2x_2& & & =& 30,\\ & 2x_0& +& 1x_1& +& 3x_2& +& 1x_3& \ge & 15,\\ & & & 2x_1& & & +& 3x_3& \le & 25,\end{array}\end{array}$$

under the bounds

$$\begin{array}{r}\begin{array}{ccccc}0& \le & x_0& \le & \mathrm{\infty },\\ 0& \le & x_1& \le & 10,\\ 0& \le & x_2& \le & \mathrm{\infty },\\ 0& \le & x_3& \le & \mathrm{\infty }.\end{array}\end{array}$$

We start our implementation in Fusion importing the relevant modules, i.e.

#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

Next we declare an optimization model creating an instance of the Model class:

  Model::t M = new Model("lo1"); auto _M = finally([&]() { M->dispose(); });

For this simple problem we are going to enter all the linear coefficients directly:

  auto A1 = new_array_ptr<double, 1>({ 3.0, 1.0, 2.0, 0.0 });
  auto A2 = new_array_ptr<double, 1>({ 2.0, 1.0, 3.0, 1.0 });
  auto A3 = new_array_ptr<double, 1>({ 0.0, 2.0, 0.0, 3.0 });
  auto c  = new_array_ptr<double, 1>({ 3.0, 1.0, 5.0, 1.0 });

The variables appearing in problem (7.1) can be declared as one $4$-dimensional variable:

  Variable::t x = M->variable("x", 4, Domain::greaterThan(0.0));

At this point we already have variables with bounds $0\le x_i\le \mathrm{\infty }$, because the domain is applied element-wise to the entries of the variable vector. Next, we impose the upper bound on $x_1$:

  M->constraint(x->index(1), Domain::lessThan(10.0));

The linear constraints can now be entered one by one using the dot product of our variable with a coefficient vector:

  M->constraint("c1", Expr::dot(A1, x), Domain::equalsTo(30.0));
  M->constraint("c2", Expr::dot(A2, x), Domain::greaterThan(15.0));
  M->constraint("c3", Expr::dot(A3, x), Domain::lessThan(25.0));

We end the definition of our optimization model setting the objective function in the same way:

  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(c, x));

Finally, we only need to call the Model.solve method:

  M->solve();

The solution values can be attained with the method Variable.level.

  auto sol = x->level();
  std::cout << "[x0,x1,x2,x3] = " << (*sol) << "\n";

Listing 7.1 Fusion implementation of model (7.1). Click here to download.

#include <iostream>
#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  auto A1 = new_array_ptr<double, 1>({ 3.0, 1.0, 2.0, 0.0 });
  auto A2 = new_array_ptr<double, 1>({ 2.0, 1.0, 3.0, 1.0 });
  auto A3 = new_array_ptr<double, 1>({ 0.0, 2.0, 0.0, 3.0 });
  auto c  = new_array_ptr<double, 1>({ 3.0, 1.0, 5.0, 1.0 });

  // Create a model with the name 'lo1'
  Model::t M = new Model("lo1"); auto _M = finally([&]() { M->dispose(); });

  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );

  // Create variable 'x' of length 4
  Variable::t x = M->variable("x", 4, Domain::greaterThan(0.0));

  // Create constraints
  M->constraint(x->index(1), Domain::lessThan(10.0));
  M->constraint("c1", Expr::dot(A1, x), Domain::equalsTo(30.0));
  M->constraint("c2", Expr::dot(A2, x), Domain::greaterThan(15.0));
  M->constraint("c3", Expr::dot(A3, x), Domain::lessThan(25.0));

  // Set the objective function to (c^t * x)
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(c, x));

  // Solve the problem
  M->solve();

  // Get the solution values
  auto sol = x->level();
  std::cout << "[x0,x1,x2,x3] = " << (*sol) << "\n";
}