CS164 KKT Conditions Theory

KKT Conditions Theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# KKT Conditions Theory\n",
    "\n",
    "1. Write out an expression for the total potential energy of the system, as a function of variables $x_1$ and $x_2$, and constants $w$ and $L$. This is the objective function that we wish to minimize.\n",
    "\n",
    "$f(x_1, x_2) = \\frac 1 2(k_1x_1^2 + k_2(x_2-x_1)^2 + k_3(L-x_2)^2)$\n",
    "\n",
    "2. Write out the three inequality constraints that prevent the masses from overlapping or from impinging on each wall.\n",
    "\n",
    "$$x1 \\geq w/2$$\n",
    "$$x2 \\geq 3w/2$$\n",
    "$$x2\\leq L - w/2$$\n",
    "\n",
    "3. Formulate the KKT conditions for the optimization problem, defining three Lagrange multipliers $λ_i$\n",
    "\n",
    "$$L(x_1, x_2, \\lambda_1, \\lambda_2, \\lambda_3) = \\frac 1 2 (k_1x_1^2 + k_2(x_2-x_1)^2 + k_3*(L-x_2)^2) + \\lambda_1(w/2 - x_1) + \\lambda_2(3w/2 - x_2) + \\lambda_3(x_2 - L + w/2)$$\n",
    "$$L_{x_1} = k_1x_1 - x_2 + 2x_1 - 1= 0$$\n",
    "$$L_{x_2} = k_2x_2 + k_3L + k_3x_2 - x_1 - \\lambda_2 + \\lambda_3= 0$$\n",
    "$$L_{\\lambda_1}= w/2 - x_1 = 0$$\n",
    "$$L_{\\lambda_2}= 3w/2 - x_2 = 0$$\n",
    "$$L_{\\lambda_3} = x_2 - L + w/2 = 0$$\n",
    "\n",
    "4. When all of the multipliers are zero, solve for the mass positions that minimize the potential energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sympy import Eq, symbols, diff, Function\n",
    "from sympy.solvers import solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, l1, l2, l3, w, k1, k2, k3, L = symbols('x1, x2, l1, l2, l3, w, k1, k2, k3, L')\n",
    "\n",
    "Lag = Function('Lag')\n",
    "\n",
    "Lag = 0.5*(k1*x1**2 + k2*(x2-x1)**2 + k3*(L-x2)**2) + l1*((w/2) - x1) + l2*(((3*w)/2) - x2) + l3*(x2 - L + (w/2))\n",
    "\n",
    "eq1 = Eq(diff(Lag, x1), 0)\n",
    "eq2 = Eq(diff(Lag, x2), 0)\n",
    "ct1 = Eq((w/2) - x1, 0)\n",
    "ct2 = Eq(((3*w)/2) - x2, 0)\n",
    "ct3 = Eq(x2 - L + (w/2), 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{x1: -(k2*(L*k3 + l2 - l3) + l1*(k2 + k3))/(k2**2 - (k1 + k2)*(k2 + k3)),\n",
       " x2: -(k2*l1 + (k1 + k2)*(L*k3 + l2 - l3))/(k2**2 - (k1 + k2)*(k2 + k3))}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2], [x1, x2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Now consider the combinations of nonzero multipliers that can be physically realized - these correspond to the active constraints. How many such combinations are there? What are the physical positions of the masses in each of these cases?\n",
    "\n",
    "There are 8 combinations of nonzero multipliers, of which at most 3 can be global minima, depending on how we set our parameters.\n",
    "\n",
    "Some of the physical positions determined by these cases correspond to the minimum energy position being when the masses are at the extremes of the segment of length $L$, or one is pushed against the other, furthest away from their side's wall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{x1: 0.5*w,\n",
       " x2: (L*k3 + 0.5*k2*w + l2 - l3)/(k2 + k3),\n",
       " l1: -(k2*(L*k3 + l2 - l3) + 0.5*w*(k2**2 - (k1 + k2)*(k2 + k3)))/(k2 + k3)}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct1], [x1, x2, l1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{x1: (1.5*k2*w + l1)/(k1 + k2),\n",
       " x2: 1.5*w,\n",
       " l2: -(k2*l1 + 1.5*w*(k2**2 - (k1 + k2)*(k2 + k3)) + (k1 + k2)*(L*k3 - l3))/(k1 + k2)}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct2], [x1, x2, l2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{x1: (0.5*k2*(2.0*L - w) + l1)/(k1 + k2),\n",
       " x2: L - 0.5*w,\n",
       " l3: (k2*l1 + 0.5*(2.0*L - w)*(k2**2 - (k1 + k2)*(k2 + k3)) + (k1 + k2)*(L*k3 + l2))/(k1 + k2)}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct3], [x1, x2, l3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{x1: 0.5*w,\n",
       " x2: 1.5*w,\n",
       " l1: 0.5*w*(k1 - 2.0*k2),\n",
       " l2: -L*k3 + k2*w + 1.5*k3*w + l3}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct1, ct2], [x1, x2, l1, l2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{x1: 0.5*w,\n",
       " x2: L - 0.5*w,\n",
       " l1: -L*k2 + 0.5*k1*w + k2*w,\n",
       " l3: -L*k2 + k2*w + 0.5*k3*w + l2}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct1, ct3], [x1, x2, l1, l3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct2, ct3], [x1, x2, l2, l3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([eq1, eq2, ct1, ct2, ct3], [x1, x2, l1, l2, l3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. In the case of nonzero multipliers, we will now try to interpret their meaning. Considering that the restoring force exerted by each spring is equal to the spring constant times its length (Hooke’s Law), give a physical meaning to the multipliers in the two equations that arise from stationarity. HINT: Newton’s second law for masses at rest states that the sum of the forces on each mass must be zero. Think of physically how you would make the masses stay on their constraint boundaries.\n",
    "\n",
    "The multipliers are forces. At these lowest-energy positions returned by the solutions to the equations derived from KKT, we can speculate that the multipliers would indicate the values of the forces acting on the system if it was left at those points. For physical systems, the Lagrangian method returns the equations of motion in a system, and so the multipliers might just be the forces acting on each of the blocks."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}