{
 "metadata": {
  "name": "",
  "signature": "sha256:78d0399fcdd9547fc8d70f4f87a4a351fdc82b8db5387566ce4a2e097d644778"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "DCP Errors with Constraints"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Created by David Stonestrom 6/4/2014.\n",
      "\n",
      "<br><br>\n",
      "The purpose of this document is to demonstrate debugging CVXPY DCP errors resulting from constraints.  \n",
      "\n",
      "<br>Let's look at an optimization problem with a simple objective and two constraints.  \n",
      "$$ \\textbf{maximize } f\\left(x\\right) = \\sum_{i=1}^{n}{x_i}$$\n",
      "\n",
      "$$\\text{subject to } Ax \\succeq b$$\n",
      "\n",
      "$$\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\sum_{i = 1}^{n}{e^{x_i}} = c$$\n",
      "\n",
      "<br>\n",
      "where $x \\in \\mathbb{R}^n$ is the variable and $A \\in \\mathbb{R}^{n \\times m}_{+}$, $b \\in \\mathbb{R}^m$, and $c \\in \\mathbb{R}$ are the problem data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cvxpy as cvx\n",
      "import numpy\n",
      "\n",
      "m = 3; # number of inequality constraints\n",
      "n = 6; # number of variables\n",
      "\n",
      "# generate some fake data\n",
      "numpy.random.seed(0) # for repeatibility\n",
      "A = numpy.random.random([m,n])\n",
      "b = numpy.random.randn(m,1) + 0.5\n",
      "c = 100 * numpy.random.random(1)\n",
      "\n",
      "#objective\n",
      "x = cvx.Variable(n,1)\n",
      "objectiveExpression = cvx.sum(x)\n",
      "\n",
      "# writing the constraints individually with names will make debugging them easier\n",
      "constraint_1 = A*x >= b\n",
      "constraint_2 = cvx.sum(cvx.exp(x)) == c\n",
      "\n",
      "# try the problem as written\n",
      "objective = cvx.Maximize(objectiveExpression)\n",
      "constraints = [constraint_1, constraint_2]\n",
      "problem = cvx.Problem(objective, constraints)\n",
      "\n",
      "# Try solving the problem. Print error if any.\n",
      "try:\n",
      "    problem.solve()\n",
      "except Exception as e:\n",
      "    print e"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Problem does not follow DCP rules.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Diagnosing the DCP error"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All right, lets look at the three components of the problem."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Objective:    \", objective.is_dcp()\n",
      "print \"Constraint 1: \", constraint_1.is_dcp()\n",
      "print \"Constraint 2: \", constraint_2.is_dcp()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Objective:     True\n",
        "Constraint 1:  True\n",
        "Constraint 2:  False\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Okay, we can dig into constraint 2 further by examining the curvature on each side."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Constraint 2 curvatures:\"\n",
      "print constraint_2.lh_exp.curvature\n",
      "print constraint_2.rh_exp.curvature"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Constraint 2 curvatures:\n",
        "CONVEX\n",
        "CONSTANT\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since the second constraint is an equality constraint, it needs each side to be AFFINE or CONSTANT in order to be DCP compliant.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Change of variables"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can make the equality constraint DCP compliant with a change of variables $y_i = e^{x_i}$.  Note that $x_i = ln\\left(y_i\\right)$.  The optimization problem is now:\n",
      "\n",
      "$$ \\textbf{maximize } g\\left(y\\right) = \\sum_{i=1}^{n}{log\\left(y\\right)}$$\n",
      "\n",
      "$$\\text{subject to } A\\log\\left(y\\right) \\succeq b$$\n",
      "\n",
      "$$\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\sum_{i = 1}^{n}{y} = c$$\n",
      "\n",
      "\n",
      "Lets look at how this transforms the constraints.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = cvx.Variable(n,1) # y = exp(x)\n",
      "constraint_1 = A * cvx.log(y) >= b\n",
      "constraint_2 = cvx.sum(y) == c\n",
      "\n",
      "print \"constraint 1:\"\n",
      "print \"left hand side:  \", constraint_1.lh_exp.curvature\n",
      "print \"right hand side: \", constraint_1.rh_exp.curvature\n",
      "print \"DCP: \", constraint_1.is_dcp()\n",
      "\n",
      "print \"\\nconstraint 2:\"\n",
      "print \"left hand side:  \", constraint_2.lh_exp.curvature\n",
      "print \"right hand side: \", constraint_2.rh_exp.curvature\n",
      "print \"DCP: \", constraint_2.is_dcp()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "constraint 1:\n",
        "left hand side:   CONSTANT\n",
        "right hand side:  UNKNOWN\n",
        "DCP:  False\n",
        "\n",
        "constraint 2:\n",
        "left hand side:   AFFINE\n",
        "right hand side:  CONSTANT\n",
        "DCP:  True\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So we've fixed the second constraint, but broken the first one.  In this case, the variables are small enough to print out the full expression for the first constraint."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Right hand side:\\n\", constraint_1.rh_exp\n",
      "print \"\\nLeft hand side:\\n\", constraint_1.lh_exp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Right hand side:\n",
        "[[ 0.5488135   0.71518937  0.60276338  0.54488318  0.4236548   0.64589411]\n",
        " [ 0.43758721  0.891773    0.96366276  0.38344152  0.79172504  0.52889492]\n",
        " [ 0.56804456  0.92559664  0.07103606  0.0871293   0.0202184   0.83261985]] * log(var3)\n",
        "\n",
        "Left hand side:\n",
        "[[ 0.94386323]\n",
        " [ 0.83367433]\n",
        " [ 1.99407907]]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that the right and left hand sides are switched from how the constraint was entered; CVXPY does this in order to represent all constraints as $\\text{LHS } \\leq \\text{ RHS}$.  \n",
      "\n",
      "The matrix multiplying $\\ln\\left(y\\right)$ is all positive, so the right hand side should be concave.  However, CVXPY requires the user to specify when a coefficient has a particular sign if that sign matters to the convexity of the expression.  This is done with parameters."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Aparam = cvx.Parameter(m,n,nonneg=True)\n",
      "Aparam.value = A\n",
      "\n",
      "constraint_1 = Aparam * cvx.log(y) >= b\n",
      "\n",
      "print \"constraint 1:\"\n",
      "print \"left hand side:  \", constraint_1.lh_exp.curvature\n",
      "print \"right hand side: \", constraint_1.rh_exp.curvature\n",
      "print \"DCP: \", constraint_1.is_dcp()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "constraint 1:\n",
        "left hand side:   CONSTANT\n",
        "right hand side:  CONCAVE\n",
        "DCP:  True\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that both constraints are good, lets check the objective and solve the problem"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "objectiveExpression = cvx.sum(cvx.log(y))\n",
      "print \"objective Curvature: \", objectiveExpression.curvature"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "objective Curvature:  CONCAVE\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Final solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All right, everything should work now."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "constraints = [constraint_1, constraint_2]\n",
      "objective = cvx.Maximize(objectiveExpression)\n",
      "problem = cvx.Problem(objective, constraints)\n",
      "\n",
      "print \"DCP rules test: \", problem.is_dcp()\n",
      "try: \n",
      "    print \"solving... \"\n",
      "    problem.solve()\n",
      "\n",
      "    # print some stuff\n",
      "    print problem.status\n",
      "    print \"Optimal value: \", problem.value\n",
      "    print \"\\nResidual of equality constraint: \", abs(sum(y.value) - c)\n",
      "    print \"\\nResidual of inequality constraints:\\n\", Aparam.value.dot(numpy.log(y.value)) - b\n",
      "    print \"\\ny:\\n\", y.value\n",
      "    print \"\\nx:\\n\", numpy.log(y.value) \n",
      "    \n",
      "except Exception as e:\n",
      "    print(e)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "DCP rules test:  True\n",
        "solving... \n",
        "optimal"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Optimal value:  4.01095525084\n",
        "\n",
        "Residual of equality constraint:  [[ 0.00134525]]\n",
        "\n",
        "Residual of inequality constraints:\n",
        "[[  1.44312355e+00]\n",
        " [  1.82027513e+00]\n",
        " [  9.80742961e-04]]\n",
        "\n",
        "y:\n",
        "[[ 2.08613192]\n",
        " [ 2.3558799 ]\n",
        " [ 1.71047705]\n",
        " [ 1.72276186]\n",
        " [ 1.66771014]\n",
        " [ 2.28582697]]\n",
        "\n",
        "x:\n",
        "[[ 0.73531159]\n",
        " [ 0.85691429]\n",
        " [ 0.53677231]\n",
        " [ 0.54392874]\n",
        " [ 0.51145151]\n",
        " [ 0.82672787]]\n"
       ]
      }
     ],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}