# 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 \leq \sum_{j=0}^{n-1} a_{kj} x_j \leq u_k^c,\quad k=0,\ldots ,m-1,$

and the bounds

$l_j^x \leq x_j \leq u_j^x, \quad j=0,\ldots ,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{split}c = \left[ \begin{array}{c} c_0 \\ \vdots \\ c_{n-1} \end{array} \right],\end{split}$
• $$c^f$$ — fixed term in the objective,

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

$\begin{split}A = \left[ \begin{array}{ccc} a_{0,0} & \cdots & a_{0,(n-1)} \\ \vdots & \cdots & \vdots \\ a_{(m-1),0} & \cdots & a_{(m-1),(n-1)} \end{array} \right],\end{split}$
• $$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{split}\begin{array} {lccccccccl} \mbox{maximize} & 3 x_0 & + & 1 x_1 & + & 5 x_2 & + & 1 x_3 & & \\ \mbox{subject to} & 3 x_0 & + & 1 x_1 & + & 2 x_2 & & & = & 30, \\ & 2 x_0 & + & 1 x_1 & + & 3 x_2 & + & 1 x_3 & \geq & 15, \\ & & & 2 x_1 & & & + & 3 x_3 & \leq & 25, \end{array}\end{split}$

under the bounds

$\begin{split}\begin{array}{ccccc} 0 & \leq & x_0 & \leq & \infty , \\ 0 & \leq & x_1 & \leq & 10, \\ 0 & \leq & x_2 & \leq & \infty ,\\ 0 & \leq & x_3 & \leq & \infty . \end{array}\end{split}$

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\leq x_i\leq \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";
}