Primal-Dual_Interior_Point_Method.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Computational_project_part2.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "uOGoiPZjzeQV"
      },
      "source": [
        "def primal_dual_interior_point(c,Q,A,b,x,pi,z,eta):\n",
        "    \"\"\" \n",
        "    Description: \n",
        "    this function uses primal-dual interior method to \n",
        "    solve convex quadratic constrained optimization\n",
        "    Input:\n",
        "    c: coefficient of linear part\n",
        "    Q: coefficient of quadratic part\n",
        "    A: coefficient of linear constraint\n",
        "    b: value of linear constraint\n",
        "    x, pi, z: column vectors of initial infeasible solution\n",
        "    eta: damping parameter\n",
        "    Output: \n",
        "    Optimal solution of x,pi,z\n",
        "    \"\"\"\n",
        "    # import library\n",
        "    import numpy as np\n",
        "    # set tolerance\n",
        "    tol=1E-8\n",
        "    # number of variables\n",
        "    var_num=len(x)\n",
        "    # calcualte stopping criteria\n",
        "    c_1=np.linalg.norm(np.dot(A,x)-b)\n",
        "    c_2=np.linalg.norm(-np.dot(Q,x)+np.dot(A.T,pi)+z-c)\n",
        "    c_3=np.dot(x.T,z)\n",
        "    # initialize counter \n",
        "    counter=0\n",
        "    # print the information of initial infeasible solution \n",
        "    print(\"=====Iteration{}=====\".format(counter))\n",
        "    print('x: {}'.format(x))\n",
        "    print('pi: {}'.format(pi))\n",
        "    print('z: {}'.format(z))\n",
        "    print('tau: {}'.format('-'))\n",
        "    print('primal problem: {}'.format(np.dot(c.T,x)+0.5*np.dot(x.T,Q).dot(x)))\n",
        "    print('dual problem: {}'.format(np.dot(b.T,pi)-0.5*np.dot(x.T,Q).dot(x)))\n",
        "    print('Residual of primal: {}'.format(c_1))\n",
        "    print('Residual of dual: {}'.format(c_2))\n",
        "    print('x.T dot x: {}'.format(c_3))\n",
        "    # check stopping criteria \n",
        "    while (c_1 > tol) | (c_2>tol) | (c_3>tol):\n",
        "        # calculate residules for the infeasible solution\n",
        "        r_p=np.dot(A,x)-b\n",
        "        r_d=-np.dot(Q,x)+np.dot(A.T,pi)+z-c\n",
        "        XZe=np.dot(np.diag(x),np.diag(z).dot(np.ones((var_num))))\n",
        "        # stack the coefficient matrix\n",
        "        coef_1=np.vstack((-Q,A,np.diag(z)))\n",
        "        coef_2=np.vstack((A.T,np.zeros((coef_1.shape[0]-var_num,A.shape[0]))))\n",
        "        coef_3=np.vstack((np.identity(var_num),np.zeros((A.shape[0],var_num)),np.diag(x)))\n",
        "        coef=np.hstack((coef_1,coef_2,coef_3))\n",
        "        # stack the residuals\n",
        "        r=np.vstack((-r_d.reshape(-1,1),-r_p.reshape(-1,1),-XZe.reshape(-1,1)))\n",
        "        # solve for affine scaling directions\n",
        "        d_aff=np.linalg.solve(coef,r)\n",
        "        # get direction for each part\n",
        "        d_x_aff=d_aff[:var_num]\n",
        "        d_pi_aff=d_aff[var_num:coef.shape[0]-var_num]\n",
        "        d_z_aff=d_aff[-var_num:]\n",
        "        # find step length of the affine scaling direction\n",
        "        # create a list to store possible selections\n",
        "        selections=[1]\n",
        "        # iterate through possible selections and append result to the list\n",
        "        for i in range(var_num):\n",
        "            if d_x_aff[i] < 0:\n",
        "                selections.append((-x[i]/d_x_aff[i])[0])\n",
        "            if d_z_aff[i] < 0:\n",
        "                selections.append((-z[i]/d_z_aff[i])[0])\n",
        "        # find the minimum value in possible selections \n",
        "        selections=np.array(selections)\n",
        "        alpha_aff=np.min(selections)\n",
        "        # calculate the duality measure \n",
        "        y=np.dot(x.T,z)/var_num\n",
        "        y_aff=np.dot((x+alpha_aff*d_x_aff.reshape(-1)).T,(z+alpha_aff*d_z_aff.reshape(-1)))/var_num\n",
        "        # calculate the centering parameter \n",
        "        tau=(y_aff/y)**3\n",
        "        # stack the adjusted residuals\n",
        "        r_l=-XZe.reshape(-1,1)+np.diag(d_x_aff.reshape(-1)).dot(np.diag(d_z_aff.reshape(-1))).dot(np.ones((var_num,1)))+tau*y*np.ones((var_num,1))\n",
        "        r_adjust=np.vstack((-r_d.reshape(-1,1),-r_p.reshape(-1,1),r_l))\n",
        "        # solve for search direction\n",
        "        d=np.linalg.solve(coef,r_adjust)\n",
        "        # get direction for each part\n",
        "        d_x=d[:var_num]\n",
        "        d_pi=d[var_num:coef.shape[0]-var_num]\n",
        "        d_z=d[-var_num:]\n",
        "        # create a list to store possible selections\n",
        "        selections=[1]\n",
        "        # iterate through possible selections and append result to the list\n",
        "        for i in range(var_num):\n",
        "            if d_x[i] < 0:\n",
        "                selections.append((-eta*x[i]/d_x[i])[0])\n",
        "            if d_z[i] < 0:\n",
        "                selections.append((-eta*z[i]/d_z[i])[0])\n",
        "        # find the minimum value in possible selections \n",
        "        selections=np.array(selections)\n",
        "        alpha=np.min(selections)\n",
        "        # calcualte new iterates\n",
        "        x=x+alpha*d_x.reshape(-1) \n",
        "        pi=pi+alpha*d_pi.reshape(-1)\n",
        "        z=z+alpha*d_z.reshape(-1)\n",
        "        # calcualte stopping criteria\n",
        "        c_1=np.linalg.norm(np.dot(A,x)-b)\n",
        "        c_2=np.linalg.norm(-np.dot(Q,x)+np.dot(A.T,pi)+z-c)\n",
        "        c_3=np.dot(x.T,z)\n",
        "        # print the information of each iteration \n",
        "        print(\"=====Iteration{}=====\".format(counter+1))\n",
        "        print('x: {}'.format(x))\n",
        "        print('pi: {}'.format(pi))\n",
        "        print('z: {}'.format(z))\n",
        "        print('tau: {}'.format(tau))\n",
        "        print('primal problem: {}'.format(np.dot(c.T,x)+0.5*np.dot(x.T,Q).dot(x)))\n",
        "        print('dual problem: {}'.format(np.dot(b.T,pi)-0.5*np.dot(x.T,Q).dot(x)))\n",
        "        print('Residual of primal: {}'.format(c_1))\n",
        "        print('Residual of dual: {}'.format(c_2))\n",
        "        print('x.T dot x: {}'.format(c_3))\n",
        "        # update counter\n",
        "        counter+=1\n",
        "    return x,pi,z\n"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YY-IP-wDe8tX"
      },
      "source": [
        "**Test the function with a quadratic constraint problem:**\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KMx7mttm4B4Y"
      },
      "source": [
        "import numpy as np\n",
        "# define the qudratic problem\n",
        "c=np.array([0,0,0])\n",
        "Q=np.array([[0.02778,0.00387,0.00021],[0.00387,0.01112,-0.0002],[0.00021,-0.0002,0.00115]])\n",
        "A=np.array([[0.1073,0.0737,0.0627],[1,1,1]])\n",
        "b=np.array([0.0650,1])\n",
        "# define the initial infeasible solution\n",
        "x=np.array([1,1,1])\n",
        "pi=np.array([1,1])\n",
        "z=np.array([1,1,1])"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "z8WLyCZA47GZ",
        "outputId": "839b5e31-5ed1-4e0c-e3b6-46a1ff6ecaa1"
      },
      "source": [
        "x_optimal,pi_optimal,z_optimal=primal_dual_interior_point(c,Q,A,b,x,pi,z,0.95)"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "=====Iteration0=====\n",
            "x: [1 1 1]\n",
            "pi: [1 1]\n",
            "z: [1 1 1]\n",
            "tau: -\n",
            "primal problem: 0.023905\n",
            "dual problem: 1.0410949999999999\n",
            "Residual of primal: 2.0079675520286675\n",
            "Residual of dual: 3.5772208141656563\n",
            "x.T dot x: 3\n",
            "=====Iteration1=====\n",
            "x: [0.05       0.551804   0.62333086]\n",
            "pi: [-6.81522183 -0.04482545]\n",
            "z: [1.01338108 0.7850784  0.70481909]\n",
            "tau: 0.005894069989513771\n",
            "primal problem: 0.0019956154242216344\n",
            "dual problem: -0.48981048232987884\n",
            "Residual of primal: 0.22603174811563292\n",
            "Residual of dual: 0.4026785558385087\n",
            "x.T dot x: 0.9232139443968299\n",
            "=====Iteration2=====\n",
            "x: [0.0025     0.50919873 0.61635703]\n",
            "pi: [-20.16128226   1.02286603]\n",
            "z: [1.27549507 0.60039681 0.3738499 ]\n",
            "tau: 0.046299944505073375\n",
            "primal problem: 0.001602622878477464\n",
            "dual problem: -0.2892199388894904\n",
            "Residual of primal: 0.12856590827897474\n",
            "Residual of dual: 0.2290418700357039\n",
            "x.T dot x: 0.5393350479333184\n",
            "=====Iteration3=====\n",
            "x: [1.25000000e-04 3.44997068e-01 7.11183476e-01]\n",
            "pi: [-45.90249572   2.90202283]\n",
            "z: [2.08323231 0.54264975 0.03485061]\n",
            "tau: 0.01629424841263047\n",
            "primal problem: 0.0009037069012638179\n",
            "dual problem: -0.08254310281537955\n",
            "Residual of primal: 0.056529852805135124\n",
            "Residual of dual: 0.1007086822055196\n",
            "x.T dot x: 0.2122581533496141\n",
            "=====Iteration4=====\n",
            "x: [5.55787214e-04 2.34345151e-01 7.76453087e-01]\n",
            "pi: [-35.13799665   2.21395942]\n",
            "z: [1.56921533 0.38985217 0.00174253]\n",
            "tau: 0.01753409849697006\n",
            "primal problem: 0.0006162051484637589\n",
            "dual problem: -0.07062657082799122\n",
            "Residual of primal: 0.011399256889168196\n",
            "Residual of dual: 0.020307927271412913\n",
            "x.T dot x: 0.09358510909277533\n",
            "=====Iteration5=====\n",
            "x: [5.24183749e-04 2.08462037e-01 7.91631453e-01]\n",
            "pi: [-9.17933441  0.57696261]\n",
            "z: [4.09608493e-01 1.02351999e-01 8.71265233e-05]\n",
            "tau: 0.0016046036806465882\n",
            "primal problem: 0.0005694677262191366\n",
            "dual problem: -0.020263590256899396\n",
            "Residual of primal: 0.0006201341514880303\n",
            "Residual of dual: 0.0011047772121793105\n",
            "x.T dot x: 0.02162018839573776\n",
            "=====Iteration6=====\n",
            "x: [5.20649138e-04 2.07057629e-01 7.92453800e-01]\n",
            "pi: [-0.4928369  0.0317996]\n",
            "z: [2.20972772e-02 6.70150534e-03 4.35632616e-06]\n",
            "tau: 1.7083048286599498e-06\n",
            "primal problem: 0.0005671542407902004\n",
            "dual problem: -0.0008019511996499861\n",
            "Residual of primal: 3.2205319511149676e-05\n",
            "Residual of dual: 5.7374203664394875e-05\n",
            "x.T dot x: 0.0014025549190781544\n",
            "=====Iteration7=====\n",
            "x: [6.00284476e-04 2.06674214e-01 7.92732595e-01]\n",
            "pi: [-0.02216286  0.00226355]\n",
            "z: [1.10486386e-03 1.51914806e-03 3.80579177e-06]\n",
            "tau: 0.006193037954916246\n",
            "primal problem: 0.000566653064275297\n",
            "dual problem: 0.00025631218477287633\n",
            "Residual of primal: 7.1219037850719995e-06\n",
            "Residual of dual: 1.2687766010303139e-05\n",
            "x.T dot x: 0.00031764893918182597\n",
            "=====Iteration8=====\n",
            "x: [0.03386353 0.07178969 0.89434678]\n",
            "pi: [-0.13819711  0.00961943]\n",
            "z: [6.61548860e-03 1.31618183e-03 6.67819303e-05]\n",
            "tau: 0.5038995274998375\n",
            "primal problem: 0.0005074275977462009\n",
            "dual problem: 0.00012918931662751429\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 3.6707339412293e-18\n",
            "x.T dot x: 0.0003782382811186882\n",
            "=====Iteration9=====\n",
            "x: [0.03342589 0.07356413 0.89300998]\n",
            "pi: [-0.02757912  0.00273556]\n",
            "z: [1.62447465e-03 6.58090914e-05 1.29174972e-05]\n",
            "tau: 0.0013175325953778468\n",
            "primal problem: 0.0005067973790760373\n",
            "dual problem: 0.00043612123124838354\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 2.6090726772098588e-18\n",
            "x.T dot x: 7.067614782765513e-05\n",
            "=====Iteration10=====\n",
            "x: [0.03079954 0.08421276 0.8849877 ]\n",
            "pi: [-0.01150857  0.00172353]\n",
            "z: [8.78698054e-04 3.29045457e-06 5.41482773e-06]\n",
            "tau: 0.021155361147097925\n",
            "primal problem: 0.0005038046592752499\n",
            "dual problem: 0.00047167200529676\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 8.481439493802682e-19\n",
            "x.T dot x: 3.213265397849015e-05\n",
            "=====Iteration11=====\n",
            "x: [0.02733624 0.09825487 0.87440888]\n",
            "pi: [0.00290613 0.00080916]\n",
            "z: [2.02287777e-04 1.64522729e-07 2.88166410e-07]\n",
            "tau: 0.011344229659299319\n",
            "primal problem: 0.0005019268608487994\n",
            "dual problem: 0.0004961289328350118\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 6.590379024195219e-19\n",
            "x.T dot x: 5.79792801378745e-06\n",
            "=====Iteration12=====\n",
            "x: [0.02646032 0.10180632 0.87173335]\n",
            "pi: [0.00658834 0.00057457]\n",
            "z: [3.06212358e-05 8.22613643e-09 2.76272305e-08]\n",
            "tau: 0.00034754163950934146\n",
            "primal problem: 0.0005018247405877753\n",
            "dual problem: 0.0005009895717107544\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 2.214368944449765e-19\n",
            "x.T dot x: 8.351688770211488e-07\n",
            "=====Iteration13=====\n",
            "x: [0.02631524 0.10239458 0.87129018]\n",
            "pi: [0.00719818 0.0005357 ]\n",
            "z: [2.20734779e-06 4.11306821e-10 1.71719962e-09]\n",
            "tau: 1.6254682940515944e-06\n",
            "primal problem: 0.0005018223551552063\n",
            "dual problem: 0.0005017627299759822\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 3.204491494080824e-19\n",
            "x.T dot x: 5.962517922405285e-08\n",
            "=====Iteration14=====\n",
            "x: [0.02630455 0.10243791 0.87125754]\n",
            "pi: [0.00724311 0.00053284]\n",
            "z: [1.14029064e-07 2.05653411e-11 8.74423751e-11]\n",
            "tau: 6.428991229583575e-10\n",
            "primal problem: 0.0005018223427293778\n",
            "dual problem: 0.0005018192649545963\n",
            "Residual of primal: 0.0\n",
            "Residual of dual: 4.824168370506119e-19\n",
            "x.T dot x: 3.077774781534848e-09\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}