```{math} \DeclareMathOperator{\sign}{sgn} \newcommand{\CO}{\textbf{\rm conv}} \newcommand{\RR}{{\mathcal R}} \newcommand{\RE}{\mathbb{R}} \newcommand{\TR}{\text{T}} \newcommand{\beq}{\begin{equation}} \newcommand{\eeq}{\end{equation}} \newcommand{\bmat}{\left[\begin{array}} \newcommand{\emat}{\end{array}\right]} \newcommand{\bsmat}{\left[\begin{smallmatrix}} \newcommand{\esmat}{\end{smallmatrix}\right]} \newcommand{\barr}{\begin{array}} \newcommand{\earr}{\end{array}} \newcommand{\bsm}{\begin{smallmatrix}} \newcommand{\esm}{\end{smallmatrix}} ``` # Quadratic program A quadratic programming (QP) optimization problem is described as: $$\begin{align} \text{minimize }&\frac{1}{2}x^\text{T}Px+q^\text{T}x+r\nonumber\\ \text{subject to }&Gx \preccurlyeq h \nonumber\\ &Ax = b \end{align}$$ $P$ is assumed to be positive semidefinite. The feasible set of QP is a polygon and the objective function is a convex quadratic function. If the objective function is quadratic and the constraints include quadratic constraints, then we have a quadratically constrained quadratic program (QCQP): $$\begin{align} \text{minimize }&\frac{1}{2}x^\text{T}P_0x+q_0^\text{T}x+r_0\nonumber\\ \text{subject to }&\frac{1}{2}x^\text{T}P_ix+q_i^\text{T}x+r_i\leq 0,\ i=1\cdots,m \nonumber\\ &Ax = b \end{align}$$ where $P_i$ for $i=0,\cdots,m$ are positive semidefinite. *Example:* Consider the set of linear equations $Ax=b$ for the case when $A$ has more rows than columns. Finding an $x$ where the equality is exactly satisfied is in general impossible. However, there is a solution for an $x$ that minimizes the cost function $e^\text{T}e$ where $e=Ax-b$. The solution is even analytic and it can be written as: \begin{equation} x^\star = (A^\text{T}A)^{-1}A^\text{T} \end{equation} However, after adding linear constraints on $x$, the optimization problem does not have an analytic solution: $$\begin{align} \text{minimize }&(Ax-b)^\text{T}(Ax-b)\nonumber\\ \text{subject to }&Gx \preccurlyeq h \nonumber\\ \end{align}$$ As a numerical example, consider: $$\begin{equation}A = \bmat{cc} 1 & 1 \\ 2 & 1\\ 3 & 2 \emat,\ b=\bmat{c} 2\\ 3 \\ 4\emat\end{equation}$$ The analytical answer to $Ax=b$ is computed as: A = np.array([[1,1],[2,1],[3,2]]) b = np.array([[2], [3], [4]]) xs = np.dot(np.linalg.pinv(A),b) print(xs) [[ 1. ] [ 0.66666667]] A similar result can be reached by solving the following QP problem: A = co.matrix(A) b = co.matrix(b) x = cp.Variable(2) I = np.identity(3) objective = cp.Minimize( cp.quad_form(A*x-b, I) ) p = cp.Problem(objective) The optimal value of the objective function is: result = p.solve() print(result) 0.333333326323 and the optimal solution is: print(x.value) [ 1.00e+00] [ 6.67e-01] Now, we can add linear constraints and find the optimal solution by solving the QP problem: x = cp.Variable(2) objective = cp.Minimize( b.T*b - 2*b.T*A*x + cp.quad_form(x, A.T*A) ) constraints = [ -0.9 <= x <= 0.9] p = cp.Problem(objective, constraints) The optimal cost function is equal to: result = p.solve() print(result) 0.33838157573 which is more than what it was without the linear constraints. The optimal solution is equal to: print(x.value) [ 9.00e-01] [ 8.18e-01] *Example (Linear Program with a Stochastic Objective Function):* Consider a random vector $c$ and the following LP problem: $$\begin{align} \text{minimize }&c^\text{T}x\nonumber\\ \text{subject to }& Gx \preccurlyeq h \nonumber\\ & Ax = b \end{align}$$ Assume that $c$ is a random vector with the normal distribution of $\mathcal{N}(\bar c,\Sigma)$. Also we assume that $x$, the unknown vector, is deterministic. With this assumptions, the objective function $c^\text{T}x$ is a normal random variable with mean ${\bar c}^\text{T}x$ and variance $x^\text{T}\Sigma x$. One way to formulate the problem so that it is practically solveable is to set the objective function as: \begin{equation} {\bar c}^\text{T}x+\gamma x^\text{T}\Sigma x \end{equation} where $\gamma\geq 0$. This objective function is called the risk-sensitive cost and $\gamma$ is call the risk-aversion parameter. The larger $\gamma$ is, the more the uncertainty of the original objective function is penalized and it thus leads to a more certain result. With this approach, the problem is formulated as the following deterministic LP: $$\begin{align} \text{minimize }&{\bar c}^\text{T}x+\gamma x^\text{T}\Sigma x\nonumber\\ \text{subject to }& Gx \preccurlyeq h \nonumber\\ & Ax = b \end{align}$$ As a numerical example, let us consider an uncertain version of ref{LP2}: $$\begin{align} \bar c=&\bmat{c} 22\\ 14.5 \emat \nonumber\\ \Sigma=&\bmat{ccc} 5 & 1\\ 1 & 4 \emat\nonumber\\ G =&\bmat{cc} -1 & -1\\ -0.5 & 1\\ 2 & -1 \emat \nonumber\\ h =&\bmat{c} 1\\2\\4\emat \end{align}$$ Now, the optimization can be solved with the following code: Sigma = co.matrix([ 5, 1, 1, 4,], (2,2)) cb = co.matrix([1, 1], (2,1)) G = co.matrix([ -1, -0.5, 2, -1, 1, -1], (3,2)) h = co.matrix([ 1, 2, 4],(3,1)) gamma = 0.5 x = cp.Variable(2) objective = cp.Minimize( cb.T * x + gamma * cp.quad_form(x, Sigma) ) constraints = [ G*x <= h ] p = cp.Problem(objective, constraints) The optimal value of the objective function is: result = p.solve() print(result) -0.184210521637 The optimal solution, in this case, is inside of the feasible set: print(x.value) [-1.58e-01] [-2.11e-01] In the following figure, the feasible set and the contours of the objective function are drawn. xp = np.array([-3, 6]) # plot the constraints plt.plot(xp, -xp-1, xp, 0.5*xp+2, xp, 2*xp-4) # Draw the feasible set (filled triangle) path = mpl.path.Path([[4, 4], [1, -2], [-2, 1], [4, 4]]) patch = mpl.patches.PathPatch(path, facecolor='green') # Add the triangle to the plot plt.gca().add_patch(patch) plt.xlabel('x') plt.ylabel('y') plt.title('Feasible Set of a Linear Program') plt.xlim(-3,6) plt.ylim(-5,5) delta = 0.025 xc = np.arange(-3.0, 6.0, delta) yc = np.arange(-5.0, 6.0, delta) X, Y = np.meshgrid(xc, yc) Z = cb[0]*X + cb[1]*Y + Sigma[0,0]*X*X + 2*Sigma[0,1]*X*Y + Sigma[1,1]*Y*Y plt.contour(X, Y, Z) X, Y = np.meshgrid(x, y) plt.show() ![png](cvxguide_73_0.png) *Example (Distance between polyhedra):* Consider the following two polyhedra: \begin{equation} \mathcal{P}_1=\{x| A_1x \preccurlyeq b_1\},\ \mathcal{P}_2=\{x| A_2x \preccurlyeq b_2\} \end{equation} The distance between $\mathcal{P}_1$ and $\mathcal{P}_2$ is defined as: \begin{equation} \text{d}(\mathcal{P}_1,\mathcal{P}_2)=\inf\{\|x_1-x_2\|_2 | x_1\in\mathcal{P}_1,\ x_2\in\mathcal{P}_2\} \end{equation} This ditance can computed using the following QP problem: $$\begin{align} \text{minimize }& \|x_1-x_2\|_2^2\nonumber\\ \text{subject to }&A_1x \preccurlyeq b_1 \nonumber\\ &A_2x \preccurlyeq b_2 \end{align}$$ As a numerical example, consider the following polygons: # Draw the triangle path = mpl.path.Path([[4, 4], [1, -2], [-2, 1], [4, 4]]) patch = mpl.patches.PathPatch(path, facecolor='green') # Add the triangle to the plot plt.gca().add_patch(patch) # Draw the rectangle path = mpl.path.Path([[-3, 4], [-2, 4], [-2, 2], [-3, 2], [-3, 4]]) patch = mpl.patches.PathPatch(path, facecolor='green') # Add the rectangle to the plot plt.gca().add_patch(patch) plt.xlabel('x') plt.ylabel('y') plt.xlim(-4,6) plt.ylim(-3,5) plt.show() ![png](cvxguide_75_0.png) The distance between these two polygons is computed with the following QP optimization problem: x1 = cp.Variable(2) x2 = cp.Variable(2) I = np.identity(2) # Triangle A1 = co.matrix([ -1, -0.5, 2, -1, 1, -1], (3,2)) b1 = co.matrix([ 1, 2, 4],(3,1)) # Rectangle A2 = co.matrix([ -1, 1, 0, 0, 0, 0, -1, 1], (4,2)) b2 = co.matrix([ 3, -2, -2, 4],(4,1)) objective = cp.Minimize( cp.quad_form(x1-x2, I) ) constraints = [ A1*x1<= b1, A2*x2<=b2] p = cp.Problem(objective, constraints) The distance between the two polygons is: result=p.solve() print(result) 0.799999994069 The correspondin point in the triangle is: print(x1.value) [-1.60e+00] [ 1.20e+00] and the corresponding point in the rectangle is: print(x2.value) [-2.00e+00] [ 2.00e+00]