snippets.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Calibration\n",
    "if mach == 0.6:\n",
    "    dt_ref = 0.0001 * subsampling\n",
    "    shed_period = 0.03\n",
    "\n",
    "elif mach == 0.1:\n",
    "    dt_ref = 0.001 * subsampling\n",
    "    shed_period = 0.18\n",
    "\n",
    "elif mach == 0.01:\n",
    "    dt_ref = 0.02 * subsampling\n",
    "    shed_period = 1.84\n",
    "\n",
    "dim_factor = shed_period * 0.16 # strouhal\n",
    "alpha = 0.2 # calibration parameter\n",
    "\n",
    "# optimal parameters \n",
    "b1_opt = b1\n",
    "b2_opt = b2\n",
    "L1_opt = L1\n",
    "L2_opt = L2\n",
    "Q_opt = Q_\n",
    "\n",
    "# computing delta denominator\n",
    "delta_denominator = 0\n",
    "delta_denominator += np.sum(np.square(b1))\n",
    "delta_denominator += np.sum(np.square(b2))\n",
    "delta_denominator += np.sum(np.square(L1))\n",
    "delta_denominator += np.sum(np.square(L2))\n",
    "delta_denominator += np.sum(np.square(Q_))\n",
    "\n",
    "# computing epsilon denominator\n",
    "epsilon_denominator = 0\n",
    "err_zero = 0\n",
    "t_cal = shed_period * 5\n",
    "steps = int(t_cal / (dt * dim_factor))\n",
    "for i in range(steps):\n",
    "    sample = int((dt * dim_factor)/dt_ref) * i\n",
    "    a = coeffs[:num,sample]\n",
    "    a_dot = np.empty_like(a)\n",
    "    for k in range(a_dot.shape[0]):\n",
    "        a_dot[k] = (nu * b1[k] + b2[k] + np.inner((nu*L1[k,:]+L2[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q_[k]),np.expand_dims(a,1)))\n",
    "    a_opt = (coeffs[:num,sample+1] - coeffs[:num,sample]) / dt_ref\n",
    "    epsilon_denominator += np.sum(np.square(a_opt - a_dot))\n",
    "    err_zero += np.sum(np.square(a_opt))\n",
    "\n",
    "def epsilon():\n",
    "\n",
    "    global Q_, L1, L2, b1, b2, nu, num\n",
    "    global Q_opt, L1_opt, L2_opt, b1_opt, b2_opt\n",
    "    \n",
    "    # calibration period \n",
    "    t_cal = shed_period * 5\n",
    "    steps = int(t_cal / (dt * dim_factor))\n",
    "    \n",
    "    numerator = 0\n",
    "    \n",
    "    for i in range(steps):\n",
    "        sample = int((dt * dim_factor)/dt_ref) * i\n",
    "\n",
    "        # optimal activation\n",
    "        a = coeffs[:num,sample]\n",
    "        # construct hypothesized activation flow\n",
    "        a_dot_opt = np.empty_like(a)\n",
    "        for k in range(a_dot.shape[0]):\n",
    "            a_dot_opt[k] = (nu * b1_opt[k] + b2_opt[k] + np.inner((nu*L1_opt[k,:]+L2_opt[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q_opt[k]),np.expand_dims(a,1)))\n",
    "        \n",
    "        # optimal activation flow\n",
    "        a_opt = (coeffs[:num,sample+1] - coeffs[:num,sample]) / dt_ref\n",
    "\n",
    "        numerator += np.sum(np.square(a_opt - a_dot_opt))\n",
    "\n",
    "    return numerator / epsilon_denominator\n",
    "\n",
    "# Delta term\n",
    "def delta():\n",
    "    # normalized distance to hypothesized model\n",
    "    global Q_, L1, L2, b1, b2, nu, num\n",
    "    global Q_opt, L1_opt, L2_opt, b1_opt, b2_opt\n",
    "\n",
    "    numerator = 0\n",
    "    denom = 0\n",
    "\n",
    "    # manual computation of components\n",
    "    numerator += np.sum(np.square(b1_opt - b1))\n",
    "    numerator += np.sum(np.square(b2_opt - b2))\n",
    "    numerator += np.sum(np.square(L1_opt - L1))\n",
    "    numerator += np.sum(np.square(L2_opt - L2))\n",
    "    for i in range(num):\n",
    "        numerator += np.sum(np.square(Q_opt[i] - Q_[i]))\n",
    "\n",
    "    return numerator / delta_denominator\n",
    "\n",
    "# gradient of f (opt_parameters)\n",
    "def grad():\n",
    "    global num\n",
    "\n",
    "    derivative = np.zeros(num)\n",
    "    for k in range(num):\n",
    "        derivative[k] += 0\n",
    "\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(8,8))\n",
    "ax.plot(coeffs[0,:],coeffs[1,:],linewidth=1, linestyle=\"dashed\", label=\"Mode 0 over Mode 1\")\n",
    "ax.plot(coeffs[2,:],coeffs[3,:],linewidth=1, linestyle=\"dashed\", label=\"Mode 2 over Mode 3\")\n",
    "#ax.plot(coeffs[4,:],coeffs[5,:],linewidth=1, linestyle=\"dashed\", label=\"Mode 4 over Mode 5\")\n",
    "ax.set_xlabel(\"Mode x\",size=15)\n",
    "ax.set_ylabel(\"Mode y\",size=15)\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(8,8))\n",
    "ax.plot(sol.y[0,-3000:len(t)],sol.y[1,-3000:len(t)],linewidth=1, linestyle=\"dashed\", label=\"Mode 0 over Mode 1\")\n",
    "ax.plot(sol.y[2,-3000:len(t)],sol.y[3,-3000:len(t)],linewidth=1, linestyle=\"dashed\", label=\"Mode 2 over Mode 3\")\n",
    "#ax.plot(sol.y[4,-1000:len(t)],sol.y[5,-1000:len(t)],linewidth=1, linestyle=\"dashed\", label=\"Mode 4 over Mode 5\")\n",
    "ax.set_xlabel(\"Mode x\",size=15)\n",
    "ax.set_ylabel(\"Mode y\",size=15)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruct flows\n",
    "rec_len = 100\n",
    "ssubsampling=20\n",
    "rec = np.repeat(np.expand_dims(q_avg,1),rec_len,1)\n",
    "for i in range(num):\n",
    "    rec += np.outer(pod_modes[:,i],sol.y[i,-100*ssubsampling::ssubsampling])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = np.zeros(200)\n",
    "for k in tqdm(range(200)):\n",
    "    nu = 3.0 + k * 0.005 \n",
    "\n",
    "    sol = solve_ivp(galerkin_system,(t0,tmax),a0,method='LSODA',t_eval=sampling_span)\n",
    "    \n",
    "    # Limit cycle -> activation limits\n",
    "    max_activ = np.empty(num)\n",
    "    min_activ = np.empty(num)\n",
    "\n",
    "    for i in range(num):\n",
    "        max_activ[i] = coeffs[i,:].max()\n",
    "        min_activ[i] = coeffs[i,:].min()\n",
    "\n",
    "    # Limit cycle -> activation limits of ROM\n",
    "    max_activ_rom = np.empty(num)\n",
    "    min_activ_rom = np.empty(num)\n",
    "\n",
    "    for i in range(num):\n",
    "        max_activ_rom[i] = sol.y[i,-int(len(t)/2):].max()\n",
    "        min_activ_rom[i] = sol.y[i,-int(len(t)/2):].min()\n",
    "\n",
    "    diff[k] = np.sum(np.divide(abs(max_activ - max_activ_rom) + abs(min_activ - min_activ_rom),1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "if controlled==True:\n",
    "    def galerkin_system(t,a):\n",
    "        # set global variables references\n",
    "        global Q_, L1, L2, b1, b2, nu\n",
    "        T = 1 / (0.165 * mach * sos) \n",
    "        T1 = T / 2 # based on doubled vortex shedding frequency f=1/T\n",
    "        T2 = T * 2 \n",
    "        A = 1\n",
    "        a_dot = np.empty_like(a)\n",
    "        gamma = np.sin((2*np.pi*t*0.005)/(T1)+A*np.pi*np.sin((2*np.pi*t*0.005)/(T2)))\n",
    "        for k in range(a_dot.shape[0]):\n",
    "            a_dot[k] = (nu * b1[k] + b2[k] + np.inner((nu*L1[k,:]+L2[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q_[k]),np.expand_dims(a,1))) + gamma * (nu * d1[k] + d2[k] + np.inner(g[k],a) + gamma * f[k]) #+ h[k] * dgamma\n",
    "        return a_dot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Limit Cycles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Limit cycle -> activation limits\n",
    "max_activ = np.empty(num)\n",
    "min_activ = np.empty(num)\n",
    "\n",
    "for i in range(num):\n",
    "    max_activ[i] = coeffs[i,:].max()\n",
    "    min_activ[i] = coeffs[i,:].min()\n",
    "\n",
    "# Limit cycle -> activation limits of ROM\n",
    "max_activ_rom = np.empty(num)\n",
    "min_activ_rom = np.empty(num)\n",
    "\n",
    "for i in range(num):\n",
    "    max_activ_rom[i] = sol.y[i,-int(len(sampling_span)/2):].max()\n",
    "    min_activ_rom[i] = sol.y[i,-int(len(sampling_span)/2):].min()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Penalty Stabilization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### bc penalty term\n",
    "# velocity on boundary \n",
    "boundary_weight = np.zeros(3*n)\n",
    "\n",
    "# compute ds\n",
    "rad = mesh.points_pol[-200,0]\n",
    "ds = 2 * np.pi * rad / 199\n",
    "\n",
    "# set boundary weights for skalar product\n",
    "boundary_weight[n-200:n] = ds\n",
    "boundary_weight[2*n-200:2*n] = ds\n",
    "boundary_weight[3*n-200:3*n] = ds * mach_weight\n",
    "\n",
    "# state based freestream condition\n",
    "q_inf = np.mean(q_fluc,1) # time averaged boundary\n",
    "\n",
    "# bc\n",
    "G1 = np.empty((num,num))\n",
    "G2 = np.empty(num)\n",
    "for i in range(num):\n",
    "    for j in range(i,num):\n",
    "        G1[i,j] = skalar_product(pod_modes[:,i],pod_modes[:,j],boundary_weight)\n",
    "        G1[j,i] = G1[i,j]\n",
    "    G2[i] = skalar_product(pod_modes[:,i],q_inf,boundary_weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    # doesnt work properly -> unphysical\n",
    "elif stabilization==\"penalty\":\n",
    "        if mach == 0.6:\n",
    "            # tbd\n",
    "            tau = 6\n",
    "        elif mach == 0.1:\n",
    "            # tbd\n",
    "            tau = 4\n",
    "        elif mach == 0.01:\n",
    "            # tbd\n",
    "            tau = 8\n",
    "\n",
    "\n",
    "    elif stabilization==\"penalty\":\n",
    "        if mach == 0.6:\n",
    "            # tbd\n",
    "            tau = 6 \n",
    "        elif mach == 0.1:\n",
    "            # tbd\n",
    "            tau = 4\n",
    "        elif mach == 0.01:\n",
    "            # tbd\n",
    "            tau = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimization Loop\n",
    "tmp = 1\n",
    "cost = 0\n",
    "counter = 0\n",
    "maxIter = 10\n",
    "alpha = 0.001\n",
    "while(abs(cost) - abs(tmp) < epsilon):\n",
    "    # max iteration loop break\n",
    "    if counter > maxIter:\n",
    "        break\n",
    "    # Setup AD Taping\n",
    "    RealReverse.setActive() \n",
    "    if counter==0:\n",
    "        design_vars = np.array([RealReverse(1.0) for _ in range(design_dim)])\n",
    "    else:\n",
    "        design_vars =  np.array([design_var for design_var in design_vars])                                              \n",
    "    # set tape variables\n",
    "    for i in range(design_dim):                                                 \n",
    "        design_vars[i].registerInput()\n",
    "\n",
    "    # copy previous cost\n",
    "    tmp = cost\n",
    "    # compute control depending on design variables\n",
    "    Nt = int(tmax/dt)\n",
    "    control, control_deriv = gamma(Nt*2,t0,tmax,dt/2,T,design_vars)\n",
    "    t = t0\n",
    "    a_tmp = a0\n",
    "    a_vec = np.zeros((num,Nt))\n",
    "    for i in range(Nt-2):\n",
    "        k1 = dt * galerkinSystem(t,a_tmp,control[2*i],control_deriv[2*i])\n",
    "        k2 = dt * galerkinSystem(t+dt/2,a_tmp+k1/2,control[2*i+1],control_deriv[2*i+1])\n",
    "        k3 = dt * galerkinSystem(t+dt/2,a_tmp+k2/2,control[2*i+1],control_deriv[2*i+1])\n",
    "        k4 = dt * galerkinSystem(t+dt,a_tmp+k3,control[2*i+2],control_deriv[2*i+2])\n",
    "        k = (k1 + 2 * k2 + 2 * k3 + k4) / 6\n",
    "        a_tmp = a_tmp + k\n",
    "        a_vec[:,i] = a_tmp\n",
    "        t = t + dt \n",
    "    # compute cost function\n",
    "    cost = costFunction(a_vec,control[::2])\n",
    "    # computing gradient with AD\n",
    "    cost.registerOutput()\n",
    "    RealReverse.setPassive()\n",
    "    cost.setGradient(1.0)\n",
    "    RealReverse.evaluate()\n",
    "    for i in range(design_dim):                                                 \n",
    "        design_vars[i] = design_vars[i] - alpha * design_vars[i].getGradient()\n",
    "\n",
    "    print(design_vars[7].getGradient())\n",
    "    # loop output\n",
    "    cost = 0\n",
    "    counter += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameterized control function\n",
    "def gamma(Nt,t0,tmax,dt,T,design_vars): \n",
    "    # empty control and control deriv arrays\n",
    "    control = [0.0 for i in range(Nt)]\n",
    "    control_deriv = [0.0 for i in range(Nt)]\n",
    "    # loop timesteps\n",
    "    for i in range(Nt):\n",
    "        t = t0 + dt * i\n",
    "        # loop design vars\n",
    "        for j in range(design_dim): \n",
    "            control[i] = control[i] + design_vars[j] * np.sin(j * np.pi * t / (T * 2))\n",
    "            control_deriv[i] = control_deriv[i] + design_vars[j] * np.cos(j * np.pi * t / (T * 2)) * j * np.pi / (T * 2)\n",
    "    \n",
    "        control[i] = control[i] + gamma_start * (1 - t / tmax)\n",
    "        control_deriv[i] = control_deriv[i] - gamma_start * (1 / tmax)\n",
    "    \n",
    "    # return control and control derivative \n",
    "    return control, control_deriv \n",
    "\n",
    "def costFunction(a,gamma):\n",
    "    cost = 0\n",
    "    for i in range(num):\n",
    "        cost = np.sum(np.sum(np.square(a),0) + np.exp(np.square(gamma) - np.square(gamma_max)) * dt)\n",
    "    return cost\n",
    "\n",
    "# Galerkin system\n",
    "def galerkinSystem(t,a,gamma,dgamma):\n",
    "    a_dot = np.empty_like(a)\n",
    "    # iterate dof\n",
    "    for k in range(a.shape[0]):\n",
    "        a_dot[k] = nu * b1[k] + b2[k] + np.inner((nu * L1[k,:]+L2[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q_[k]),np.expand_dims(a,1)) + gamma * (nu * d1[k] + d2[k] + np.inner(g[k],a) + gamma * f[k]) + h[k] * dgamma\n",
    "    # return activation derivative\n",
    "    return a_dot\n",
    "\n",
    "def solver(f,control_function,design_vars,interval,t0,y0,T,dt=dt):\n",
    "    t0 = interval[0]\n",
    "    tmax = interval[1]\n",
    "    Nt = int(tmax/dt)\n",
    "    control, control_deriv = control_function(Nt*2,t0,tmax,dt/2,T,design_vars)\n",
    "    t = t0\n",
    "    y_vec = np.zeros((num,Nt))\n",
    "    y_vec[:,0] = y0\n",
    "    for i in range(Nt-1):\n",
    "        y = y_vec[:,i]\n",
    "        k1 = dt * f(t,y,control[2*i],control_deriv[2*i])\n",
    "        k2 = dt * f(t+dt/2,y+k1/2,control[2*i+1],control_deriv[2*i+1])\n",
    "        k3 = dt * f(t+dt/2,y+k2/2,control[2*i+1],control_deriv[2*i+1])\n",
    "        k4 = dt * f(t+dt,y+k3,control[2*i+2],control_deriv[2*i+2])\n",
    "        k = (k1+2*k2+2*k3+k4)/6\n",
    "        y_vec[:,i+1] = y_vec[:,i] + k\n",
    "        t = t + dt  \n",
    "    return y_vec\n",
    "\n",
    "def costateProblem(t,phi,a,gamma):\n",
    "    phi_dot = np.zeros_like(phi)\n",
    "    for k in range(phi.shape[0]):\n",
    "        A = np.zeros(a.shape[0])\n",
    "        for i in range(a.shape[0]):\n",
    "            tmp = 0\n",
    "            for j in range(a.shape[0]):\n",
    "                tmp += (Q_[k][i,j] + Q_[k][j,i]) * a[j]\n",
    "            A[i] = - (nu * L1[k][i] + L2[k][i] + tmp + g[k,i] * gamma)\n",
    "        phi_dot[k] = np.inner(A,phi) - 0.5 * a[k]\n",
    "    return phi_dot\n",
    "\n",
    "def coStateSolver(f,control_function,design_vars,interval,t0,y0,T,dt=dt):\n",
    "    t0 = interval[0]\n",
    "    tmax = interval[1]\n",
    "    Nt = int(tmax/dt)\n",
    "    control, control_deriv = control_function(Nt*2,t0,tmax,dt/2,T,design_vars)\n",
    "    t = t0\n",
    "    y_vec = np.zeros((num,Nt))\n",
    "    y_vec[:,0] = y0\n",
    "    for i in range(Nt-1):\n",
    "        y = y_vec[:,i]\n",
    "        k1 = dt * f(t,y,control[2*i],control_deriv[2*i])\n",
    "        k2 = dt * f(t+dt/2,y+k1/2,control[2*i+1],control_deriv[2*i+1])\n",
    "        k3 = dt * f(t+dt/2,y+k2/2,control[2*i+1],control_deriv[2*i+1])\n",
    "        k4 = dt * f(t+dt,y+k3,control[2*i+2],control_deriv[2*i+2])\n",
    "        k = (k1+2*k2+2*k3+k4)/6\n",
    "        y_vec[:,i+1] = y_vec[:,i] + k\n",
    "        t = t + dt  \n",
    "    return y_vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'num' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-328dd30f759a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;31m# LIST OF LOCAL VARIABLES\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0mdof\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnum\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m \u001b[0mpoly_len\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdof\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdof\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'num' is not defined"
     ]
    }
   ],
   "source": [
    "# single dof calibration\n",
    "\n",
    "# LIST OF LOCAL VARIABLES\n",
    "\n",
    "dof = num\n",
    "poly_len = int((dof+1)*(dof+2)/2) \n",
    "\n",
    "\n",
    "# construct polinomial\n",
    "polynomial = np.ones((dof,poly_len))  \n",
    "# mapping polynomials to galerkin coefficients\n",
    "for k in range(dof):\n",
    "    polynomial[k][0] = nu * b1[k] + b2[k]\n",
    "    polynomial[k][1:dof+1] = nu * L1[k,:] + L2[k,:]\n",
    "    tmp = 0\n",
    "    for i in range(dof):\n",
    "        tmp2 = 0\n",
    "        for j in range(dof-1,i-1,-1):\n",
    "            if i==j:\n",
    "                polynomial[k][dof+1+tmp+tmp2] = Q_[k][i,i]\n",
    "                tmp2 += 1\n",
    "            else:\n",
    "                polynomial[k][dof+1+tmp+tmp2] = Q_[k][i,j] + Q_[k][j,i]\n",
    "                tmp2 += 1\n",
    "        tmp += dof - i \n",
    "\n",
    "orig_polynomial = polynomial "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if stabilization==\"calibration\":\n",
    "    alpha = 0.1 # calibration parameter\n",
    "\n",
    "class galerkin_calibration:\n",
    "    \n",
    "    def __init__(self, b1, b2, L1, L2, Q_, activations, T, dim):\n",
    "        \n",
    "        # fixed parameters\n",
    "        self.alpha = 0.2\n",
    "        self.nu = 1.716e-5\n",
    "        \n",
    "        # optimization parameters\n",
    "        self.dof = dim\n",
    "        tmp = self.dof*(self.dof+1)*(self.dof+2)/2 \n",
    "        self.polynomial = np.ones((self.dof,tmp))    \n",
    "        \n",
    "        # mapping polynomials to galerkin coefficients\n",
    "        for i in range(self.dof):\n",
    "            self.polynomial[i][0] = nu * b1[i] + b2[i]\n",
    "            self.polynomial[i][1:self.dof+1] = nu * L1[i,:] + L2[i,:]\n",
    "            for k in range(self.dof):\n",
    "                for j in range(self.dof-1):\n",
    "                    if k==i:\n",
    "                        self.polynomial[i][self.dof+1+(k*self.dof - k)+j] = Q_[j,j]\n",
    "                    else:\n",
    "                        self.polynomial[i][self.dof+1+(k*self.dof - k)+j] = Q_[k,j] + Q_[j,k]\n",
    "        self.orig_polynomial = self.polynomial \n",
    "\n",
    "        # constructing activation mapping \n",
    "        self.activations = activations\n",
    "        \n",
    "        self.dim_factor = shed_period * 0.16        \n",
    "        self.error = np.zeros(self.dof)\n",
    "        self.distance = np.zeros(self.dof)\n",
    "        self.cost = np.zeros(self.dof)\n",
    "\n",
    "        # temporal setup        \n",
    "        self.dt = dt\n",
    "        self.calibration_time = shed_period * 10\n",
    "        self.temporal_dim_factor = shed_period * 0.16\n",
    "        self.discrete_steps = int(self.calibration_time / (self.dt * self.temporal_dim_factor))\n",
    "\n",
    "        # compute denominators of delta and epsilon\n",
    "        self.delta_denominator = np.sum(np.square(self.orig_polynomial),1)\n",
    "        self.epsilon_denominator = 0\n",
    "        \n",
    "        for t in range(self.discrete_steps):\n",
    "            ref_sample_time = int(dt/ dt_ref) * t \n",
    "            a = self.activations[:self.dof,ref_sample_time]\n",
    "            a_dot = np.empty_like(a)\n",
    "            for k in range(self.dof):\n",
    "                a_dot[k] = (nu * b1[k] + b2[k] + np.inner((nu*L1[k,:]+L2[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q_[k]),np.expand_dims(a,1)))\n",
    "            a_opt = (self.activations[:self.dof,ref_sample_time+1] - self.activations[:self.dof, ref_sample_time]) / dt_ref\n",
    "            self.epsilon_denominator += np.square(a_opt - a_dot)\n",
    "\n",
    "    def compute_error():\n",
    "        self.error = np.zeros(self.dof)\n",
    "        for t in range(self.discrete_steps):\n",
    "            ref_sample_time = int(dt/ dt_ref) * t\n",
    "            a = self.activations[:self.dof,ref_sample_time]\n",
    "            a_dot = np.empty_like(a)\n",
    "            for k in range(a_dot.shape[0]):\n",
    "                a_dot[k] = self.polynomial[k,0] + np.inner(self.polynomial[k,1:self.dof+1],a)\n",
    "                for i in range(a_dot.shape[0]):\n",
    "                    for j in range(a_dot.shape[0]-i):\n",
    "                        a_dot[k] += a[i] * a[j] * self.polynomial[i][self.dof+1+(i*self.dof - i)+j]\n",
    "            a_opt = (self.activations[:self.dof, ref_sample_time+1] - self.activations[:self.dof, ref_sample_time]) / dt_ref\n",
    "            self.error += np.square(a_opt - a_dot)\n",
    "\n",
    "    def compute_distance():\n",
    "        self.distance = np.sum(np.square(self.polynomial - self.orig_polynomial),1) \n",
    "\n",
    "    def compute_cost():\n",
    "        self.cost = (1 - self.alpha) * (self.error / self.epsilon_denominator) + self.alpha * (self.distance / self.delta_denominator)\n",
    "\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# control function\n",
    "T = shed_period\n",
    "def control(t):\n",
    "    tmp = control0 * (1 - t / T)\n",
    "    for m in range(designDim):\n",
    "        tmp += designVars[m] * np.sin((m+1) * np.pi * t * 0.5 / T)\n",
    "    return tmp\n",
    "\n",
    "def controlDerivative(t):\n",
    "    tmp = - control0 * (1 / T)\n",
    "    for m in range(designDim):\n",
    "        tmp += designVars[m] * np.cos((m+1) * np.pi * t * 0.5 / T) * m * np.pi * 0.5 / T\n",
    "    return tmp\n",
    "\n",
    "def control_galerkin_system(t, a):\n",
    "    n = a.shape[0]\n",
    "    a_dot = np.zeros(n)\n",
    "    gamma = control(t)\n",
    "    dgamma = controlDerivative(t)\n",
    "    # iterate dof\n",
    "    for k in range(n):\n",
    "        a_dot[k] = nu * b1[k] + b2[k] + np.inner((nu * L1[k,:]+L2[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q[k]),np.expand_dims(a,1)) + gamma * (nu * d1[k] + d2[k] + np.inner(g[k],a) + gamma * f[k]) + h[k] * dgamma\n",
    "    # return activation derivative\n",
    "    return a_dot\n",
    "\n",
    "sol = solve_ivp(control_galerkin_system,(t0,tmax),a0, method='LSODA',t_eval=temporalSpan)\n",
    "sol.message\n",
    "\n",
    "plt.plot(temporalSpan,[control(t) for t in temporalSpan])\n",
    "for i in range(4):\n",
    "    plt.plot(sol.y[i,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# old optimization \n",
    "\n",
    "class controlOptimization:\n",
    "\n",
    "    def __init__(self, designDim, timeInterval, dt, T, a0, phi0, control0, maxControl):\n",
    "        \n",
    "        self.designDim = designDim\n",
    "        self.designVariables = np.zeros(designDim)\n",
    "\n",
    "        # temporal setup\n",
    "        self.dt = dt\n",
    "        self.t0 = timeInterval[0]\n",
    "        self.tmax = np.ceil(timeInterval[1] / dt) * dt\n",
    "        self.T = T\n",
    "        self.Nt = int(self.tmax/dt) + 1\n",
    "\n",
    "        # initial conditions\n",
    "        self.a0 = a0\n",
    "        self.phi0 = phi0\n",
    "        self.control0 = control0\n",
    "\n",
    "        # misc\n",
    "        self.dim = len(a0)\n",
    "        self.maxControl = maxControl\n",
    "        self.N = 3 # control weight in cost \n",
    "        self.cost = 0\n",
    "\n",
    "        # optimizaiton control\n",
    "        self.epsilon = 1e-6\n",
    "        self.alpha = 1e-6\n",
    "        self.maxIter = 1000\n",
    "\n",
    "        # initialize arrays \n",
    "        self.a = np.zeros((self.dim, self.Nt))   # system trajectory\n",
    "        self.phi = np.zeros((self.dim, self.Nt)) # costate function\n",
    "        self.gamma = np.zeros(self.Nt)           # control trajectory\n",
    "        self.grad = np.zeros(self.designDim)     # gradient\n",
    "\n",
    "    # printing statements\n",
    "    def status(self,print_string):\n",
    "        print(f'{print_string:{\"-\"}<80}')\n",
    "\n",
    "    # functions for single timestep output\n",
    "    def control(self,t):\n",
    "        tmp = self.control0 * (1 - t / self.T)\n",
    "        for m in range(self.designDim):\n",
    "            tmp += self.designVariables[m] * np.sin((m+1) * np.pi * t * 0.5 / self.T)\n",
    "        return tmp\n",
    "\n",
    "    def controlDerivative(self,t):\n",
    "        tmp = - self.control0 * (1 / self.T)\n",
    "        for m in range(self.designDim):\n",
    "            tmp += self.designVariables[m] * np.cos((m+1) * np.pi * t * 0.5 / self.T) * (m+1) * np.pi * 0.5 / self.T\n",
    "        return tmp\n",
    "\n",
    "    def activation(self,t):\n",
    "        # interpolation of state trajectory a\n",
    "        if t % self.dt != 0:\n",
    "            return 0.5 * (self.a[:,int(np.ceil(t/self.dt))] + self.a[:,int(np.floor(t/self.dt))])\n",
    "        else:\n",
    "            return self.a[:,int(t/self.dt)]     \n",
    "\n",
    "    # RK45 ODE solver for forward/reverse time solve\n",
    "    def solver(self,f,interval,y0,dt):\n",
    "        # RK45 Solver -> trajectory of y\n",
    "        t = interval[0]\n",
    "        tmax = interval[1]\n",
    "        Nt = np.abs(int((tmax-t)/dt)) + 1\n",
    "        # reverse time alternative\n",
    "        alt = False\n",
    "        if (tmax < t):\n",
    "            dt *= -1\n",
    "            alt = True\n",
    "        # initial condition\n",
    "        y = np.zeros((y0.shape[0],Nt))\n",
    "        y[:,0] = y0\n",
    "        # stepping loop\n",
    "        for i in range(Nt-1):\n",
    "            k1 = dt * f(t,y[:,i])\n",
    "            k2 = dt * f(t+dt/2,y[:,i]+k1/2)\n",
    "            k3 = dt * f(t+dt/2,y[:,i]+k2/2)\n",
    "            k4 = dt * f(t+dt,y[:,i]+k3)\n",
    "            k = (k1+2*k2+2*k3+k4)/6\n",
    "            y[:,i+1] = y[:,i] + k\n",
    "            t = t + dt  \n",
    "        if alt:\n",
    "            return np.flip(y,1)\n",
    "        else:\n",
    "            return y\n",
    "    \n",
    "    # ODE system dynamics\n",
    "    def costateProblem(self, t, phi):\n",
    "        n = phi.shape[0]\n",
    "        phi_dot = np.zeros(n)\n",
    "        # time step activations\n",
    "        a = self.activation(t)\n",
    "        gamma = self.control(t)\n",
    "        # looping dimensions\n",
    "        for k in range(n):\n",
    "            # costate linear vector\n",
    "            A = np.zeros(n)\n",
    "            # computing costate vector\n",
    "            for i in range(n):\n",
    "                tmp = 0\n",
    "                for j in range(n):\n",
    "                    tmp += (Q[i][k,j] + Q[i][j,k]) * a[j]\n",
    "                A[i] = - (nu * L1[i][k] + L2[i][k] + tmp + g[i,k] * gamma)\n",
    "            phi_dot[k] = np.inner(A,phi) - 0.5 * a[k]\n",
    "        return phi_dot\n",
    "\n",
    "    def galerkinSystem(self, t, a):\n",
    "        n = a.shape[0]\n",
    "        a_dot = np.zeros(n)\n",
    "        gamma = self.control(t)\n",
    "        dgamma = self.controlDerivative(t)\n",
    "        # iterate dof\n",
    "        for k in range(n):\n",
    "            a_dot[k] = nu * b1[k] + b2[k] + np.inner((nu * L1[k,:]+L2[k,:]),a) + np.matmul(np.matmul(np.expand_dims(a,1).T,Q[k]),np.expand_dims(a,1)) + gamma * (nu * d1[k] + d2[k] + np.inner(g[k],a) + gamma * f[k]) + h[k] * dgamma\n",
    "        # return activation derivative\n",
    "        return a_dot\n",
    "\n",
    "    # discrete gradient computation   \n",
    "    def gradientFunction(self,phi):\n",
    "        gradient = np.zeros(self.designDim)\n",
    "        # gradient computation loop\n",
    "        for m in range(self.designDim):\n",
    "            for n in range(self.Nt):    \n",
    "                t = self.t0 + n * self.dt\n",
    "                gamma = self.control(t)\n",
    "                tmpTime = 0\n",
    "                for k in range(self.dim):\n",
    "                    # computing vector based A\n",
    "                    A = np.zeros(self.dim)\n",
    "                    for i in range(self.dim):\n",
    "                        tmp = 0\n",
    "                        for j in range(self.dim):\n",
    "                            tmp += (Q[i][j,k] + Q[i][k,j]) * self.a[j,n]\n",
    "                        A[i] = - (nu * L1[i][k] + L2[i][k] + tmp + g[i][k] * gamma)\n",
    "                    # computing temporal gradient step\n",
    "                    tmp = 0\n",
    "                    for i in range(self.dim):\n",
    "                        tmp  += h[k] * A[i] * phi[i,n] + g[k][i] * self.a[i,n] * phi[k,n]\n",
    "                    tmpTime += tmp + (nu * d1[k] + d2[k] + 2 * gamma * f[k]) * phi[k][n] - h[k] * 0.5 * self.a[k][n]\n",
    "                gradient[m] += self.dt * (tmpTime  + np.exp(self.N * (np.power(gamma,2) - np.power(gamma_max,2))) * 2 * self.N * gamma) * np.sin((m+1) * np.pi * t * 0.5 / self.T)\n",
    "        return gradient\n",
    "\n",
    "        # discrete gradient computation   \n",
    "    def gradientFunction(self,phi):\n",
    "        gradient = np.zeros(self.designDim)\n",
    "        # gradient computation loop\n",
    "        for m in range(self.designDim):\n",
    "            for n in range(self.Nt):    \n",
    "                t = self.t0 + n * self.dt\n",
    "                gamma = self.control(t)\n",
    "                tmpTime = 0\n",
    "                for k in range(self.dim):\n",
    "                    # computing vector based A\n",
    "                    A = np.zeros(self.dim)\n",
    "                    for i in range(self.dim):\n",
    "                        tmp = 0\n",
    "                        for j in range(self.dim):\n",
    "                            tmp += (Q[i][j,k] + Q[i][k,j]) * self.a[j,n]\n",
    "                        A[i] = - (nu * L1[i][k] + L2[i][k] + tmp + g[i][k] * gamma)\n",
    "                    # computing temporal gradient step\n",
    "                    tmp = 0\n",
    "                    for i in range(self.dim):\n",
    "                        tmp  += h[k] * A[i] * phi[i,n] + g[k][i] * self.a[i,n] * phi[k,n]\n",
    "                    tmpTime += tmp + (nu * d1[k] + d2[k] + 2 * gamma * f[k]) * phi[k][n] - h[k] * 0.5 * self.a[k][n]\n",
    "                gradient[m] += self.dt * (tmpTime  + np.exp(self.N * (np.power(gamma,2) - np.power(gamma_max,2))) * 2 * self.N * gamma) * np.sin((m+1) * np.pi * t * 0.5 / self.T)\n",
    "        return gradient\n",
    "\n",
    "    def costFunction(self):\n",
    "        # cost function\n",
    "        cost = 0\n",
    "        for n in range(self.Nt - 1):\n",
    "            t = self.t0 + n * self.dt\n",
    "            cost += self.dt * (np.sum(np.square(self.a[:,n])) + np.exp(self.N * (np.power(self.control(t),2)-np.power(self.maxControl,2)))) \n",
    "        return cost\n",
    "\n",
    "    # optimization functions\n",
    "    def costFunctionDerivatives(self):\n",
    "        costActivationDerivative = np.zeros(self.dim)\n",
    "        costControlDerivative = 0\n",
    "        for n in range(self.dim):\n",
    "            costActivationDerivative += self.dt * 0.5 * a[:,n]\n",
    "            costControlDerivative += self.dt * np.exp(self.N * (np.square(self.gamma[n]) - np.square(self.maxControl))) * self.N * 2 * self.gamma[n]\n",
    "        return [costActivationDerivative, costControlDerivative]\n",
    "\n",
    "    def sideConstraints(self):\n",
    "        constraintCost = np.zeros(self.dim)\n",
    "        for n in range(self.Nt):\n",
    "            t = self.t0 + n * self.dt\n",
    "            aDot = self.galerkinSystem(t,self.a[:,n])\n",
    "            for k in range(self.dim):\n",
    "                constraintCost[k] += self.dt * self.phi[k,n] * (aDot[k] - nu * b1[k] - b2[k] - np.inner(nu * L1[k,:] + L2[k,:] ,self.a[:,n]) - np.matmul(np.matmul(self.a[:,n],Q[k]),self.a[:,n]) - self.gamma[n] * (nu * d1[k] + d2[k]) - np.square(self.gamma[n]) - self.gamma[n] * np.inner(g[k],self.a[:,n]) - h[k] * self.controlDerivative(t))  \n",
    "        return constraintCost\n",
    "\n",
    "    def step(self):\n",
    "        # construct control trajectory\n",
    "        for n in range(self.Nt):\n",
    "            self.gamma[n] = self.control(self.t0 + n * self.dt)\n",
    "        # solving galerkin system\n",
    "        self.a = self.solver(self.galerkinSystem,[self.t0,self.tmax],self.a0,self.dt)\n",
    "        # solving costate ode dynamics in reverse time\n",
    "        self.phi = self.solver(self.costateProblem,[self.tmax,self.t0],self.phi0,self.dt)\n",
    "        # computing gradient step function\n",
    "        self.grad = self.gradientFunction(self.phi)\n",
    "        # perform update step\n",
    "        self.designVariables -= self.alpha * self.grad\n",
    "        # computing cost function\n",
    "        self.cost = self.costFunction()\n",
    "\n",
    "    def optimize(self):\n",
    "        self.status(\"Optimization started with parameters \")\n",
    "        self.status(\"Learning Rate: \" + str(self.alpha))\n",
    "        self.status(\"Iteration Limit: \" + str(self.maxIter))\n",
    "        # optimization loop\n",
    "        for s in range(self.maxIter):\n",
    "            # previouse cost\n",
    "            prevCost = self.cost\n",
    "            # computation step\n",
    "            self.step()\n",
    "            # step status\n",
    "            sys.stdout.write(\"\\r\"f'{\"Iteration: \" + str(s):{\"\"}<20}' + f'{\"Cost: \" + str(self.cost):{\"\"}<60}')\n",
    "            sys.stdout.flush()\n",
    "            # discontinuation criterion\n",
    "            if np.isnan(self.cost):\n",
    "                self.status(\"Optimization failed and stopped! \")\n",
    "                break\n",
    "            if np.abs(prevCost - self.cost) < self.epsilon:\n",
    "                self.status(\"\")\n",
    "                self.status(\"Optimization converged and finished! \")\n",
    "                break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_rk45_solver():\n",
    "    def solver(f,interval,y0,dt):\n",
    "        # RK45 Solver -> trajectory of y\n",
    "        t = interval[0]\n",
    "        tmax = interval[1]\n",
    "        Nt = np.abs(int((tmax-t)/dt)) + 1\n",
    "        # reverse time alternative\n",
    "        alt = False\n",
    "        if (tmax < t):\n",
    "            dt *= -1\n",
    "            alt = True\n",
    "        # initial condition\n",
    "        y = np.zeros((y0.shape[0],Nt))\n",
    "        y[:,0] = y0\n",
    "        # stepping loop\n",
    "        for i in range(Nt-1):\n",
    "            k1 = np.abs(dt) * f(t,y[:,i])\n",
    "            k2 = np.abs(dt) * f(t+dt/2,y[:,i]+  k1/2)\n",
    "            k3 = np.abs(dt) * f(t+dt/2,y[:,i] + k2/2)\n",
    "            k4 = np.abs(dt) * f(t+dt,y[:,i] + k3)\n",
    "            k = (k1+2*k2+2*k3+k4)/6\n",
    "            y[:,i+1] = y[:,i] + k\n",
    "            t = t + dt  \n",
    "        if alt:\n",
    "            return np.flip(y,1)\n",
    "        else:\n",
    "            return y\n",
    "\n",
    "    # comparison of custom rk45\n",
    "    #fw and bw scipy\n",
    "    solFwScipy = solve_ivp(galerkin_system,(t0,tmax),a0, method='RK45',t_eval=temporalSpan)\n",
    "    solBwScipy = solve_ivp(galerkin_system,(tmax,t0),a0, method='RK45',t_eval=np.flip(temporalSpan))\n",
    "    #fw and bw custom\n",
    "    solFwOwn = solver(galerkin_system,(t0,tmax),a0,0.0005)\n",
    "    solBwOwn = solver(galerkin_system,(tmax,t0),a0,0.0005)\n",
    "\n",
    "    fig,ax = plt.subplots(1,2,figsize=(20,6))\n",
    "    fig.tight_layout(pad=3.0)\n",
    "    for i in range(4):\n",
    "        ax[0].plot(solFwScipy.y[i,:],label=\"Mode \"+str(i))\n",
    "        ax[0].plot(solFwOwn[i,:],label=\"Mode Own \"+str(i))\n",
    "        ax[0].legend(title=\"Activation\")\n",
    "        ax[0].set_xlabel(\"time\")\n",
    "        ax[0].set_ylabel(\"amplitude\")\n",
    "\n",
    "    for i in range(4):\n",
    "        ax[1].plot(solBwScipy.y[i,:],label=\"Mode \"+str(i))\n",
    "        ax[1].plot(solBwOwn[i,:],label=\"Mode Own \"+str(i))\n",
    "        ax[1].legend(title=\"Activation\")\n",
    "        ax[1].set_xlabel(\"time\")\n",
    "        ax[1].set_ylabel(\"amplitude\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use of poincare bendixson -> outer and inner radius tending to limit cycle\n",
    "discretization = 32\n",
    "\n",
    "# specify outer ellipsoid radii\n",
    "r1 = 2*np.sqrt(np.sum(np.square(a0[:2])))\n",
    "r2 = 2*np.sqrt(np.sum(np.square(a0[2:]))) \n",
    "\n",
    "# compute angular spacing\n",
    "angles = np.linspace(0,2*np.pi,discretization)\n",
    "\n",
    "# coords in first dims\n",
    "dim1Values = np.sin(angles) * r1\n",
    "dim2Values = np.cos(angles) * r1\n",
    "dim3Values = np.sin(angles) * r2\n",
    "dim4Values = np.cos(angles) * r2\n",
    "\n",
    "# double loop to iterate all coord combinations\n",
    "for i in range(discretization):\n",
    "    for j in range(discretization):\n",
    "        \n",
    "        # construct state vector\n",
    "        state = np.array([dim1Values[i],dim2Values[i],dim3Values[j],dim4Values[j]])\n",
    "\n",
    "        # compute next state with euler fw step\n",
    "        nextState = galerkin_system(0,state) * dt + state\n",
    "\n",
    "        # compute radius changes\n",
    "        r1diff = np.sqrt(np.sum(np.square(nextState[:2]))) - r1\n",
    "        r2diff = np.sqrt(np.sum(np.square(nextState[2:]))) - r2\n",
    "\n",
    "        print(r1diff,r2diff)\n",
    "        # check if radius decreased\n",
    "        if r1diff > 0 or r2diff > 0:\n",
    "            print(\"No limit cycle in state ellipsoid!\")\n",
    "            break\n",
    "    else:\n",
    "        continue\n",
    "    break\n",
    "# equilibrium as inner radius -> instable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_optimization_results():\n",
    "    # initialize canvas\n",
    "    fig,ax = plt.subplots(3,1,figsize=(15,10))\n",
    "    ax[0].set_title(\"\")\n",
    "\n",
    "    # plot curves\n",
    "    for i in range(optimizer.dim):\n",
    "        ax[0].plot(optimizer.a[i,:],label=\"ROM Mode \"+str(i+1),c=colormap[i+1])\n",
    "        ax[0].plot(stateTrajectory[i,:],c=colormap[i+1],linestyle='dotted')\n",
    "        ax[1].plot(optimizer.phi[i,:],label=\"Costate Function \"+str(i+1),c=colormap[i+1])\n",
    "\n",
    "    for i in range(optimizer.dim):\n",
    "        ax[0].plot(stateTrajectory[i,:],c=colormap[i+1],linestyle='dashed',label=\"Reference Mode \"+str(i+1))\n",
    "    # plot control function\n",
    "    ax[2].plot(optimizer.gamma, label=\"Control Function\")\n",
    "\n",
    "    # force legends\n",
    "    ax[0].legend()\n",
    "    ax[1].legend()\n",
    "    ax[2].legend();\n",
    "\n",
    "if visual:\n",
    "    plot_optimization_results()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "5c2a26602e28ad08aeb3db5218c8af55fd3775ea6ffbecb1728e15703c31277b"
  },
  "kernelspec": {
   "display_name": "Python 3.9.5 64-bit",
   "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.9.5"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}