# Quadratic programming examples

A quadratic programming problem can be expressed in the following form:

$$
\begin{align}
\text{minimize} & \quad \frac{1}{2}x^T Q x + p^T x \\
\text{subject to} & \quad
\begin{cases}
G x \le h \\
A x = b
\end{cases}
\end{align}
$$

where:
- $p$ is a real-valued, $n$-dimensional vector;
- $Q$ is an $n \times n$-dimensional real symmetric matrix;
- $G$ is an $m \times n$-dimensional real matrix;
- $h$ is an $m$-dimensional real vector;
- $A$ is a $p \times n$-dimensional real matrix;
- $b$ is a $p$-dimensional real vector.

We will use the [qpsolvers](https://github.com/stephane-caron/qpsolvers) library, together with [numpy](https://numpy.org/), to formulate and solve a few examples of quadratic programming problems. Both libraries can be installed with `pip` as explained in the README file; once installed, we can import them as follows:


```python
import numpy as np

from qpsolvers import solve_qp
```

## Example 1 

This example is taken from [CVXOPT](https://cvxopt.org/examples/tutorial/qp) and is formulated as follows:

$$
\begin{align}
\text{minimize} & \quad f(x) = 2x_1^2 + x_2^2 + x_1x_2 + x_1 + x_2 \\
\text{subject to} & \quad
\begin{cases}
x_1 \ge 0 \\
x_2 \ge 0 \\
x_1 + x_2 = 1
\end{cases}
\end{align}
$$

We need to extract the $Q$ matrix first. We have 2 variables, so for the quadratic term we have:

$
\begin{align}
q(x) & = \dfrac{1}{2}x^T Q x \\
& = \dfrac{1}{2} \begin{bmatrix}x_1 & x_2\end{bmatrix}\begin{bmatrix}a_{11} & a_{12} \\ a_{21} & a_{22}\end{bmatrix}\begin{bmatrix}x_1 \\ x_2\end{bmatrix} \\
& = \dfrac{1}{2} \sum\limits_{i=1}^2 \sum\limits_{j=1}^2 a_{ij} x_i x_j \\
& = \dfrac{1}{2} (a_{11} x_1 x_1 + a_{12} x_1 x_2 + a_{21} x_2 x_1 + a_{22} x_2 x_2) \\
& = \dfrac{1}{2} (a_{11} x_1^2 + a_{22} x_2^2 + (a_{12} + a_{21}) x_1 x_2)
\end{align}
$

Since $Q$ is assumed to be symmetric we have that $a_{12} = a_{21}$, so that:

$
q(x) = \dfrac{1}{2} a_{11} x_1^2 + \dfrac{1}{2} a_{22} x_2^2 + a_{12} x_1 x_2
$

Now we can find the values for the elements $a_{ij}$:

$
\begin{cases}
\dfrac{1}{2} a_{11} & = 2 \implies a_{11} = 4 \\
\dfrac{1}{2} a_{22} & = 1 \implies a_{22} = 2 \\
a_{12} & = 1
\end{cases}
$

and write the $Q$ matrix as:

$Q = \begin{bmatrix}4 & 1 \\ 1 & 2\end{bmatrix}$

or, in order to show the quadratic coefficients more clearly, as:

$Q = 2 * \begin{bmatrix}2 & \dfrac{1}{2} \\ \dfrac{1}{2} & 1\end{bmatrix}$


```python
Q = 2 * np.array([[2, 0.5], [0.5, 1]])
```

The linear term is straightforward:

$
p(x) = p^T x = \begin{bmatrix}x_1 & x_2\end{bmatrix}\begin{bmatrix}a_1 \\ a_2\end{bmatrix} = a_1 x_1 + a_2 x_2
$

resulting in:

$
\begin{cases}
a_1 = 1 \\
a_2 = 1
\end{cases}
$

from which the $p$ vector can be written as:

$
p = \begin{bmatrix}1 & 1\end{bmatrix}
$


```python
p = np.array([1.0, 1.0])
```

We have two inequality constraints, $x_1 \ge 0$ and $x_2 \ge 0$. We need to write them in the form $G x \le h$, so we have:

$
\begin{align}
\begin{cases}
-x_1 &+ 0x_2 & \le 0 \\
0x_2 &- x_2 & \le 0
\end{cases}
\end{align}
$

from which we derive:

$
\begin{bmatrix}-1 & 0 \\ 0 & -1\end{bmatrix} \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} \le \begin{bmatrix}0 \\ 0 \end{bmatrix} \implies 
G = \begin{bmatrix}-1 & 0 \\ 0 & -1\end{bmatrix} , h = \begin{bmatrix}0 \\ 0 \end{bmatrix}
$


```python
G = np.array([[-1.0, 0.0], [0.0, -1.0]])
h = np.array([0.0, 0.0])
```

Finally, we have an equality constraint $x_1 + x_2 = 1$, that we can represent as:

$
\begin{bmatrix}1 & 1\end{bmatrix} \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix}1\end{bmatrix} \implies 
A = \begin{bmatrix}1 & 1\end{bmatrix} , b = \begin{bmatrix}1\end{bmatrix}
$


```python
A = np.array([1.0, 1.0])
b = np.array([1.0])
```

We now have everything we need to run the solver, which will return the optimal solution $x^* = \begin{bmatrix}x^*_1 & x^*_2\end{bmatrix}$.


```python
x_star = solve_qp(Q, p, G, h, A, b)
print(x_star)
```

    [0.25 0.75]


The `solve_qp` method uses the [quadprog](https://pypi.python.org/pypi/quadprog/) solver by default, but it can be used with other solvers as well; take a look at the [project's page](https://github.com/stephane-caron/qpsolvers#solvers) to find the list of supported solvers.

## Example 2 

This example is taken from the [Northwestern University website](https://optimization.mccormick.northwestern.edu/index.php/Quadratic_programming#Numerical_example) and is formulated as follows:

$$
\begin{align}
\text{minimize} & \quad f(x) = 3x_1^2 + x_2^2 + 2x_1x_2 + x_1 + 6x_2 + 2 \\
\text{subject to} & \quad
\begin{cases}
2x_1 + 3x_2 \ge 4 \\
x_1 \ge 0 \\
x_2 \ge 0
\end{cases}
\end{align}
$$

The constant term $2$ in the objective function can be discarded. Since there are no equality constraints, we will only have $G$ and $h$; let's keep in mind that the inequality contraints have to be expressed as $Gx \le h$.


```python
Q = 2 * np.array([[3.0, 1.0], [1.0, 1.0]])
p = np.array([1.0, 6.0])
G = np.array([[-2.0, -3.0], [-1.0, 0.0], [0.0, -1.0]])
h = np.array([-4.0, 0.0, 0.0])

x_star = solve_qp(Q, p, G, h)
print(x_star)
```

    [0.5 1. ]


Since the objective function can be expressed as $f(x) = \frac{1}{2}x^T Q x + p^T x$, its value in $x^*$ can be found as follows:


```python
f = 0.5 * x_star.T @ Q @ x_star + p.T @ x_star + 2
print(f)
```

    11.25


## Example 3

This example is taken from the [Matlab website](https://uk.mathworks.com/help/optim/ug/quadprog.html) and is formulated as follows:

$$
\begin{align}
\text{minimize} & \quad f(x) = \dfrac{1}{2} x_1^2 + x_2^2 − x_1 x_2 − 2 x_1 − 6 x_2 \\
\text{subject to} & \quad
\begin{cases}
x_1 + x_2 \le 2 \\
-x_1 + 2 x_2 \le 2 \\
2 x_1 + x_2 \le 3
\end{cases}
\end{align}
$$

This time, instead of using float numbers, we will declare the type of each variable explicitly with `dtype`:


```python
Q = np.array([[1, -1], [-1, 2]], dtype=np.float)
p = np.array([-2, -6], dtype=np.float)

G = np.array([[1, 1], [-1, 2], [2, 1]], dtype=np.float)
h = np.array([2, 2, 3], dtype=np.float)

x_star = solve_qp(Q, p, G, h)
print(x_star)
```

    [0.66666667 1.33333333]


Value of the objective function in $x^*$:


```python
f = 0.5 * x_star.T @ Q @ x_star + p.T @ x_star
print(f)
```

    -8.222222222222218


9eb557f2add12e4728a361dedc2096c9949de778
ipynb
Jupyter Notebook
qp.ipynb
nvitucci/notebook_qp
1c8489a597a1caff80c4f179328f14d3ac85de9f
MIT
qp.ipynb
nvitucci/notebook_qp
1c8489a597a1caff80c4f179328f14d3ac85de9f
MIT
qp.ipynb
nvitucci/notebook_qp
1c8489a597a1caff80c4f179328f14d3ac85de9f
MIT
Qwen/Qwen-72B
1. YES
2. YES