{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spring Optimization\n",
    "\n",
    "\n",
    "From: Sandgren, E. (1990). Nonlinear integer and discrete programming in mechanical design optimization. Transactions\n",
    "of the ASME, Journal of Mechanical Design, 112(2):223–229, June 1990. ISSN 0738-0666\n",
    "\n",
    "**Global Optimum** (according to Lampinen): 2.65856\n",
    "\n",
    "Note: In the runs below, you should be able to yield the exact global minimum, but in general, because a stochastic optimization method is being used, you will usually get close to, but not exactly to, the global minimum.  You may want to play with the control parameters of the optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from math import pi\n",
    "import PyCEGO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quasi-optimal Values of the objective function according to the literature\n",
      "Sandgren: 221838.49974522655 (typo seems likely in Lampinen)\n",
      "Chen: 2.670860274197602\n",
      "Wu and Chao: 2.6680583380916043\n",
      "Lampinen: 2.682999597110495\n"
     ]
    }
   ],
   "source": [
    "# Problem originally from Sandgren\n",
    "allowable_wire_diameters= [\n",
    "    0.009, 0.0095, 0.0104, 0.0118, 0.0128, 0.0132,\n",
    "    0.014, 0.015, 0.0162, 0.0173, 0.018, 0.020,\n",
    "    0.023, 0.025, 0.028, 0.032, 0.035, 0.041,\n",
    "    0.047, 0.054, 0.063, 0.072, 0.080, 0.092,\n",
    "    0.105, 0.120, 0.135, 0.148, 0.162, 0.177,\n",
    "    0.192, 0.207, 0.225, 0.244, 0.263, 0.283,\n",
    "    0.307, 0.331, 0.362, 0.394, 0.4375, 0.500\n",
    "]\n",
    "\n",
    "# Constants of the problem\n",
    "F_max = 1000.0 # [lb] Maximum working load\n",
    "S = 189000.0 # [psi] Allowable maximum shear stress\n",
    "l_max = 14.0 # [in] Maximum free length\n",
    "d_min = 0.2 # [in] Minimum wire diameter\n",
    "D_max = 3.0 # [in] Maximum outer diameter of the spring\n",
    "F_p = 300.0 # [lb] Preload compression force\n",
    "sigma_pm = 6.0 # [in] Maximum deflection under preload\n",
    "sigma_w = 1.25 # [in] Deflection from preload to full load\n",
    "G = 11.5e6 # [] Shear modulus of the material\n",
    "\n",
    "def CEGO_obj(x):\n",
    "    return obj([x[0].as_int(),\n",
    "                x[1].as_double(),\n",
    "                allowable_wire_diameters[x[2].as_int()]])\n",
    "\n",
    "def obj(x):\n",
    "    # Unpack the inputs\n",
    "    N,D,d = x\n",
    "    \n",
    "    # Intermediate terms\n",
    "    C_f = (4*(D/d)-1)/(4*(D/d)-4)+0.615*d/D\n",
    "    K = G*d**4/(8*N*D**3)\n",
    "    sigma_p = F_p/K\n",
    "    l_f = F_max/K + 1.05*(N+2)*d\n",
    "    \n",
    "    g = np.zeros((9,))\n",
    "    g[1] = 8*C_f*F_max*D/(pi*d**3) - S\n",
    "    g[2] = l_f - l_max \n",
    "    # g3, g4, and g5 are handled explicitly via bounds; will always be met\n",
    "    g[3] = d_min - d # will always be met\n",
    "    g[4] = D-D_max # will always be met\n",
    "    g[5] = 3.0-D/d # will always be met\n",
    "    g[6] = sigma_p-sigma_pm\n",
    "    g[7] = sigma_p+(F_max-F_p)/K + 1.05*(N+2)*d - l_f\n",
    "    g[8] = sigma_w - (F_max - F_p)/K\n",
    "    \n",
    "    s = np.ones_like(g)\n",
    "    s[7] = 1e5\n",
    "    s[8] = 1e5\n",
    "    \n",
    "    c = np.array([1 if g_i < 0 else 1+s_i*g_i for s_i,g_i in zip(s,g)])\n",
    "    assert((c >= 1).all())\n",
    "    \n",
    "    f = pi**2*D*d**2*(N+2)/4\n",
    "\n",
    "    return f*np.product(c**3)\n",
    "\n",
    "print('Quasi-optimal Values of the objective function according to the literature')\n",
    "print('Sandgren:',obj([10, 1.180701, 0.283]), \"(typo seems likely in Lampinen)\") # Sandgren (see Lampinen)\n",
    "print('Chen:',obj([9, 1.2287, 0.283])) # Chen and Tsao (see Lampinen)\n",
    "print('Wu and Chao:',obj([9, 1.227411, 0.283])) # Wu and Chao (see Lampinen)\n",
    "print('Lampinen:',obj([9, 1.2230410, 0.283])) # Lampinen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i: 0 best: 5.00286 c: 12, 1.536645, 36,  queue: 0\n",
      "i: 50 best: 3.02387 c: 18, 0.885894, 34,  queue: 0\n",
      "i: 100 best: 2.66119 c: 9, 1.224252, 35,  queue: 0\n",
      "i: 150 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 200 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 250 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 300 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 350 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 400 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 450 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 500 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 550 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 600 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 650 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 700 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 750 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 800 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 850 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 900 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 950 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 0 best: 8.68977 c: 36, 0.813757, 35,  queue: 0\n",
      "i: 50 best: 2.83485 c: 10, 1.195462, 35,  queue: 0\n",
      "i: 100 best: 2.65856 c: 9, 1.223043, 35,  queue: 0\n",
      "i: 150 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 200 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 250 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 300 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 350 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 400 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 450 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 500 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 550 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 600 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 650 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 700 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 750 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 800 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 850 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 900 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 950 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 0 best: 7.98714 c: 17, 1.555042, 37,  queue: 0\n",
      "i: 50 best: 2.70755 c: 5, 1.663264, 36,  queue: 0\n",
      "i: 100 best: 2.65963 c: 9, 1.223534, 35,  queue: 0\n",
      "i: 150 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 200 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 250 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 300 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 350 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 400 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 450 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 500 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 550 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 600 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 650 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 700 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 750 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 800 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 850 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 900 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 950 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 0 best: 3.34426 c: 4, 2.061832, 37,  queue: 0\n",
      "i: 50 best: 2.82774 c: 10, 1.192467, 35,  queue: 0\n",
      "i: 100 best: 2.65857 c: 9, 1.223044, 35,  queue: 0\n",
      "i: 150 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 200 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 250 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 300 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 350 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 400 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 450 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 500 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 550 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 600 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 650 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 700 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 750 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 800 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 850 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 900 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 950 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 0 best: 6.72479 c: 26, 1.032771, 36,  queue: 0\n",
      "i: 50 best: 2.67581 c: 9, 1.230977, 35,  queue: 0\n",
      "i: 100 best: 2.65858 c: 9, 1.223049, 35,  queue: 0\n",
      "i: 150 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 200 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 250 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 300 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 350 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 400 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 450 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 500 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 550 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 600 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 650 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 700 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 750 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 800 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 850 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 900 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n",
      "i: 950 best: 2.65856 c: 9, 1.223041, 35,  queue: 0\n"
     ]
    }
   ],
   "source": [
    "D = 3\n",
    "Nwired = len(allowable_wire_diameters)\n",
    "CEGO_bounds = [PyCEGO.Bound(1, int(l_max/d_min)), # N\n",
    "               PyCEGO.Bound(3*d_min, D_max), # D\n",
    "               PyCEGO.Bound(27, Nwired-1) # d (indices thereof)\n",
    "               ]\n",
    "\n",
    "for ocounter in range(5):\n",
    "    layers = PyCEGO.NumberishLayers(CEGO_obj, D, D*20, 1, 3)\n",
    "    layers.set_bounds(CEGO_bounds)\n",
    "    layers.set_builtin_evolver(PyCEGO.BuiltinEvolvers.differential_evolution)\n",
    "\n",
    "    objs = []\n",
    "    for counter in range(1000):\n",
    "        layers.do_generation()\n",
    "        objective, coeffs = layers.get_best()\n",
    "        if counter % 50 == 0:\n",
    "            print(layers.print_diagnostics())\n",
    "        objs.append(objective)"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}