appendix-proxqp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0fc42397-e1f9-425c-a20c-e5bf7d03f9a7",
   "metadata": {},
   "source": [
    "# The quadradic-program solver ProxQP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f9d2e01-ba9f-43f9-be98-ff97d2f210b1",
   "metadata": {},
   "source": [
    "This appendix is a very brief introduction to the solver ProxQP, dedicated to problems with quadric cost and linear (inequality and equality) constraints. We will need numpy and the ProxSuite package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d01be3d8-939b-466f-9fd1-4b4f430a3d51",
   "metadata": {},
   "outputs": [],
   "source": [
    "import proxsuite\n",
    "import numpy as np\n",
    "from numpy.linalg import eig\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4473ffeb-2bac-4469-8c0c-ee3834e60484",
   "metadata": {},
   "source": [
    "## Problem definition \n",
    "We consider the following linearly constrained quadratic program:\n",
    "$$\\min_x \\frac{1}{2} x^T H x - g^T x$$\n",
    "Under the constraints:\n",
    "$$\\underline{b} \\le x \\le \\overline{b}$$\n",
    "$$A_e x = b_e$$\n",
    "$$\\underline{b}_i \\le A_i x \\le \\overline{b}_i$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bef6de4e-8122-4984-9b8d-eeae0716596b",
   "metadata": {},
   "source": [
    "where $H$ is the cost Hessian, $g$ the cost gradient, $\\underline{b}$,$\\overline{b}$ the bounds (box constraints), $A_e$,$b_e$ define the equality constraints and $A_i,\\underline{b}_i,\\overline{b}_i$ the inequality constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "524baace-b670-402d-a502-af1c228d5e68",
   "metadata": {},
   "source": [
    "## Random construction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "119b03af-68cf-441d-9983-101cec1293c9",
   "metadata": {},
   "source": [
    "Here is a piece of code to define a problem, with the following sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ed70c5c-b839-4465-9275-41a114e88d01",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp4/generated/example_proxqp_param\n",
    "# ### TEST PARAMETERS\n",
    "NX = 20 # x dimension (search space)\n",
    "NEQ = 5 # number of equalities\n",
    "NINEQ = 3 # Number of inequalities\n",
    "WITH_BOUNDS = True # Additional bounds on x\n",
    "VERBOSE = False # Do you want to see the result?\n",
    "ACCURACY = 1e-6 # Threshold for solver stoping criteria and posterior checks\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1031e666-e184-43ad-b8aa-b1a8240fb4ef",
   "metadata": {},
   "source": [
    "Sampling is relatively straightforward. You just have to pay attention to the Hessian that should be symmetrical positive definite, and the bounds that should be sorted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5664e5a-b50d-40f2-bd85-c3f0e1fbd448",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp4/generated/example_proxqp_matrices\n",
    "# Cost\n",
    "H = np.random.rand(NX, NX) * 2 - 1\n",
    "H = H @ H.T  ### Make it positive symmetric\n",
    "g = np.random.rand(NX)\n",
    "\n",
    "Ae = np.random.rand(NEQ, NX)*2-1\n",
    "be = np.random.rand(NEQ)*2-1\n",
    "\n",
    "Ai = np.random.rand(NINEQ, NX)*2-1\n",
    "bi = np.sort(np.random.rand(NINEQ,2)*2-1,1)\n",
    "for i in range(NINEQ):\n",
    "    # Half inequalities are double bounds\n",
    "    # One quarter are pure lower\n",
    "    # One quarter are pure upper\n",
    "    r = random.randint(0,3)\n",
    "    if r==0: bi[i,0] = -1e20\n",
    "    elif r==1: bi[i,1] = 1e20\n",
    "\n",
    "bounds = np.sort(np.random.rand(NX,2)*2-1,1)+[-1,1]\n",
    "for i in range(NX):\n",
    "    # Half inequalities are double bounds\n",
    "    # One quarter are pure lower\n",
    "    # One quarter are pure upper\n",
    "    r = random.randint(0,3)\n",
    "    if r==0: bounds[i,0] = -1e20\n",
    "    elif r==1: bounds[i,1] = 1e20\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13238437-2821-442d-946e-8b67d380ee3d",
   "metadata": {},
   "source": [
    "## Solver call"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f69a64a1-c048-4c2e-bcbe-42a650a4bb67",
   "metadata": {},
   "source": [
    "The solver must first be initialized to the matrix sizes, then call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fba22f6c-d320-4b67-b2ec-3cd1886b60db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp4/generated/example_proxqp_solve\n",
    "#[x, cost, _, niter, lag, iact] =\n",
    "qp = proxsuite.proxqp.dense.QP(NX, NEQ, NINEQ, WITH_BOUNDS)\n",
    "qp.settings.eps_abs = ACCURACY/1e3\n",
    "qp.init(H, g, Ae, be, Ai, bi[:,0], bi[:,1],\n",
    "        bounds[:,0] if WITH_BOUNDS else None,\n",
    "        bounds[:,1] if WITH_BOUNDS else None)\n",
    "qp.solve()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1e0cf3a-cab5-4437-829f-569a1a9fd317",
   "metadata": {},
   "source": [
    "The results are exported through the *results* field. The solver provides the primal variable $x^*$, and the multipliers of the equality $y^*$ and inequality $z^*$ constraints (bounds and inequalities are put together in $z^*$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0993f88c-f88c-49e9-b857-7b2b04a8178f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp4/generated/example_proxqp_result\n",
    "x,y,z = qp.results.x,qp.results.y,qp.results.z\n",
    "if WITH_BOUNDS:\n",
    "    w = z[NINEQ:] # bounds\n",
    "    z = z[:NINEQ] # general inequalities\n",
    "cost = qp.results.info.objValue\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e6047b6-4e21-4e05-8c2c-5a18548eb491",
   "metadata": {},
   "source": [
    "Let's print the result to finish!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a07ca53-3a7b-4510-a8af-025269275ff6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load tp4/generated/example_proxqp_print\n",
    "print(\"Primal optimum x: {}\".format(x))\n",
    "print(\"Dual optimum (equalities) y: {}\".format(y))\n",
    "print(\"Dual optimum (ineq) z: {}\".format(z))\n",
    "print(\"Dual optimum (bounds) w: {}\".format(w))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}