diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
index 2f82207e..a3131a84 100644
--- a/examples/tutorial_1_WaveBot.ipynb
+++ b/examples/tutorial_1_WaveBot.ipynb
@@ -26,7 +26,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -35,7 +35,7 @@
     "import matplotlib.pyplot as plt\n",
     "from scipy.optimize import brute\n",
     "\n",
-    "import wecopttool as wot"
+    "import wecopttool as wot\n"
    ]
   },
   {
@@ -83,7 +83,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -92,7 +92,7 @@
     "mesh = wb.mesh(mesh_size_factor)\n",
     "fb = cpy.FloatingBody.from_meshio(mesh, name=\"WaveBot\")\n",
     "fb.add_translation_dof(name=\"Heave\")\n",
-    "ndof = fb.nb_dofs"
+    "ndof = fb.nb_dofs\n"
    ]
   },
   {
@@ -107,12 +107,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
-    "fb.show_matplotlib()\n",
-    "_ = wb.plot_cross_section(show=True)  # specific to WaveBot"
+    "# fb.show_matplotlib()\n",
+    "# _ = wb.plot_cross_section(show=True)  # specific to WaveBot\n"
    ]
   },
   {
@@ -127,22 +127,33 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Warning: Minimum wavelength in frequency spectrum (0.6939155518806638) is smaller than the minimum computable wavelength (0.8056649972775113).\n"
+     ]
+    }
+   ],
    "source": [
-    "f1 = 0.05\n",
-    "nfreq = 50\n",
-    "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n",
+    "df = 0.15\n",
+    "f1 = df\n",
+    "ifreq_start = 2\n",
+    "ifreq_end = 10 # 50\n",
+    "\n",
+    "freq = wot.frequency(df, ifreq_end)[ifreq_start:]\n",
     "\n",
     "min_computable_wavelength = fb.minimal_computable_wavelength\n",
     "g = 9.81\n",
-    "min_period = 1/(f1*nfreq)\n",
+    "min_period = 1/(df*ifreq_end)\n",
     "min_wavelength = (g*(min_period)**2)/(2*np.pi)\n",
     "\n",
     "if min_wavelength < min_computable_wavelength:\n",
     "    print(f'Warning: Minimum wavelength in frequency spectrum ({min_wavelength}) is smaller'\n",
-    "         f' than the minimum computable wavelength ({min_computable_wavelength}).')"
+    "         f' than the minimum computable wavelength ({min_computable_wavelength}).')\n"
    ]
   },
   {
@@ -157,11 +168,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "WARNING:wecopttool.core:Using the geometric centroid as the center of gravity (COG).\n",
+      "WARNING:wecopttool.core:Using the center of gravity (COG) as the rotation center for hydrostatics.\n",
+      "WARNING:capytaine.bem.problems_and_results:Mesh resolution for DiffractionProblem(body=WaveBot_immersed, omega=9.425, water_depth=inf, wave_direction=0.000, rho=1025.0):\n",
+      "The resolution of the mesh 'WaveBot' of the body 'WaveBot_immersed' might be insufficient for the wavelength λ=6.94e-01.\n",
+      "This warning appears because the largest panel of this mesh has radius 1.01e-01 > λ/8.\n",
+      "WARNING:capytaine.bem.problems_and_results:Mesh resolution for RadiationProblem(body=WaveBot_immersed, omega=9.425, water_depth=inf, radiating_dof='Heave', rho=1025.0):\n",
+      "The resolution of the mesh 'WaveBot' of the body 'WaveBot_immersed' might be insufficient for the wavelength λ=6.94e-01.\n",
+      "This warning appears because the largest panel of this mesh has radius 1.01e-01 > λ/8.\n"
+     ]
+    }
+   ],
    "source": [
-    "bem_data = wot.run_bem(fb, freq)"
+    "bem_data = wot.run_bem(fb, freq)\n"
    ]
   },
   {
@@ -181,7 +207,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -190,7 +216,7 @@
     "controller = None\n",
     "loss = None\n",
     "pto_impedance = None\n",
-    "pto = wot.pto.PTO(ndof, kinematics, controller, pto_impedance, loss, name)"
+    "pto = wot.pto.PTO(ndof, kinematics, controller, pto_impedance, loss, name)\n"
    ]
   },
   {
@@ -202,25 +228,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# PTO dynamics forcing function\n",
-    "f_add = {'PTO': pto.force_on_wec}\n",
-    "\n",
-    "# Constraint\n",
-    "f_max = 2000.0\n",
-    "nsubsteps = 4\n",
-    "\n",
-    "def const_f_pto(wec, x_wec, x_opt, waves): # Format for scipy.optimize.minimize\n",
-    "    f = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)\n",
-    "    return f_max - np.abs(f.flatten())\n",
-    "\n",
-    "ineq_cons = {'type': 'ineq',\n",
-    "             'fun': const_f_pto,\n",
-    "             }\n",
-    "constraints = [ineq_cons]"
+    "# DEBUG\n",
+    "constraints = []\n"
    ]
   },
   {
@@ -233,7 +246,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -242,7 +255,7 @@
     "    constraints=constraints,\n",
     "    friction=None,\n",
     "    f_add=f_add,\n",
-    ")"
+    ")\n"
    ]
   },
   {
@@ -267,15 +280,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [],
    "source": [
-    "amplitude = 0.0625  \n",
+    "f1 = df\n",
+    "amplitude = 0 # 0.0625\n",
     "wavefreq = 0.3\n",
     "phase = 30\n",
     "wavedir = 0\n",
-    "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
+    "waves = wot.waves.regular_wave(f1, ifreq_end, wavefreq, amplitude, phase, wavedir, ifreq_start)\n"
    ]
   },
   {
@@ -290,12 +304,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
     "obj_fun = pto.mechanical_average_power\n",
-    "nstate_opt = 2*nfreq"
+    "nstate_opt = wec.ncomponents\n"
    ]
   },
   {
@@ -315,193 +329,237 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
-    "options = {'maxiter': 200}\n",
-    "scale_x_wec = 1e1\n",
-    "scale_x_opt = 1e-3\n",
-    "scale_obj = 1e-2\n",
+    "options = {'maxiter': 100}\n",
+    "# scale_x_wec = 1.0 # 1e1\n",
+    "# scale_x_opt = 1.0 # 1e-3\n",
+    "# scale_obj = 1.0 # 1e-2\n",
+    "\n",
     "\n",
-    "results = wec.solve(\n",
-    "    waves, \n",
-    "    obj_fun, \n",
+    "scale_x_wec = 1.0\n",
+    "scale_x_opt = 1.0\n",
+    "scale_obj = 1.0\n",
+    "\n",
+    "# results  = wec.solve(\n",
+    "problem  = wec.solve(\n",
+    "    waves,\n",
+    "    obj_fun,\n",
     "    nstate_opt,\n",
-    "    optim_options=options, \n",
+    "    # x_wec_0 = x_wec*1.1,\n",
+    "    # x_opt_0 = x_opt*0.9,\n",
+    "    optim_options=options,\n",
     "    scale_x_wec=scale_x_wec,\n",
     "    scale_x_opt=scale_x_opt,\n",
     "    scale_obj=scale_obj,\n",
+    "    # use_grad = False,\n",
     "    )\n",
     "\n",
-    "opt_mechanical_average_power = results.fun\n",
-    "print(f'Optimal average mechanical power: {opt_mechanical_average_power} W')"
+    "# opt_mechanical_average_power = results.fun\n",
+    "# print(f'Optimal average mechanical power: {opt_mechanical_average_power} W')\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Analyzing results\n",
-    "We will use two post-processing functions to obtain frequency- and time-domain results for the WEC and PTO responses. The pseudospectral method gives continuous in time results. To get smoother looking plots, we specify the number of subpoints betweeen co-location points. In this case we will use 5. "
+    "## DEBUG"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "nsubsteps = 5\n",
-    "pto_fdom, pto_tdom = pto.post_process(wec, results, waves, nsubsteps=nsubsteps)\n",
-    "wec_fdom, wec_tdom = wec.post_process(results, waves, nsubsteps=nsubsteps)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
+   "execution_count": 13,
    "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.11.4\n",
+      "Optimization terminated successfully    (Exit mode 0)\n",
+      "            Current function value: 0.0004998016765365283\n",
+      "            Iterations: 8\n",
+      "            Function evaluations: 8\n",
+      "            Gradient evaluations: 8\n"
+     ]
+    }
+   ],
    "source": [
-    "The `pto.post_process` function returns `xarray.Dataset`s, which have built-in integration with PyPlot for smart plotting that automagically sets titles and formatting. We will plot the mechanical power (`mech_power`), position (`pos`), and the PTO force (`force`)."
+    "## SOLVE!!!! ##\n",
+    "from scipy.optimize import minimize\n",
+    "import scipy; print(scipy.__version__)\n",
+    "\n",
+    "optim_res = minimize(**problem)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.figure()\n",
-    "pto_tdom['power'].loc['mech',:,:].plot()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
+   "execution_count": 14,
    "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "dict_keys(['fun', 'x0', 'method', 'constraints', 'options', 'bounds', 'callback', 'jac'])"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
-    "We could similarly plot any time- or frequency-domain repsonse of the WEC or PTO. For instance, here is the PTO heave motion."
+    "problem.keys()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 15,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "'SLSQP'"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
-    "plt.figure()\n",
-    "wec_tdom['pos'].plot()"
+    "problem['method']\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 16,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'maxiter': 100, 'disp': True}"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
-    "plt.figure()\n",
-    "pto_tdom['force'].plot()"
+    "problem['options']"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": 17,
    "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
-    "Note that there are other dynamic responses available in the post-processed WEC and PTO variables (`wec_tdom`, `pto_tdom`, `wec_fdom`, `pto_fdom`). For example, the time domain PTO variable contains the following response:"
+    "problem['bounds'] == None"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [],
    "source": [
-    "pto_tdom"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## 2. Optimal control for maximum electrical power\n",
+    "# SETUP\n",
+    "# idx = 3\n",
+    "# x_wec = np.zeros(wec.ncomponents); x_wec[idx*2] = 1; x_wec[idx*2+1] = 1; x_wec[0] = 1\n",
+    "# x_opt = np.zeros(wec.ncomponents); x_opt[idx*2] = 1; x_opt[idx*2+1] = 1; x_opt[0] = 1\n",
     "\n",
-    "The rest of this tutorial will focus on optimizing for electrical power (new objective function) rather than mechanical, as this is a form of power that is usable and transportable.\n",
+    "# x = np.concatenate([x_wec, x_opt])\n",
     "\n",
-    "Since we're still dealing with the same WaveBot as in part 1, we can reuse the BEM and wave data from before. Look back at part 1 if you need a refresher on how to create these data.\n",
-    "\n",
-    "<div>\n",
-    "<img src=\"https://live.staticflickr.com/65535/52435033525_b8efc28d16_k.jpg\" width=\"1000\">\n",
-    "</div>\n",
-    "The WEC-PTO kinematics remain the same as well (unity). The major difference now is that we consider the dynamics of PTO, since they impact the electrical power and we shall not assume a lossless PTO.\n",
-    "\n",
-    "We will express the PTO's dynamics in form of a 2-port impedance model, to incoporate the dynamics of the drive-train and the dynamics of the generator.\n",
-    "The additional mechanical energy storage through the drive-train is modelled using Newton's second law and we assume a linear generator using a power-invariant park transform.\n",
-    "\n",
-    "The PTO impedance matrix components are then obtained under open-circuit  conditions, i.e.  no  load  current  or  no  WEC  velocity,  respectively."
+    "# t = wec.time"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
-    "## PTO impedance definition\n",
-    "omega = bem_data.omega.values\n",
-    "gear_ratio = 12.0\n",
-    "torque_constant = 6.7\n",
-    "winding_resistance = 0.5\n",
-    "winding_inductance = 0.0\n",
-    "drivetrain_inertia = 2.0\n",
-    "drivetrain_friction = 1.0\n",
-    "drivetrain_stiffness = 0.0\n",
-    "\n",
-    "drivetrain_impedance = (1j*omega*drivetrain_inertia + \n",
-    "                        drivetrain_friction + \n",
-    "                        1/(1j*omega)*drivetrain_stiffness) \n",
-    "\n",
-    "winding_impedance = winding_resistance + 1j*omega*winding_inductance\n",
-    "\n",
-    "\n",
-    "pto_impedance_11 = -1* gear_ratio**2 * drivetrain_impedance\n",
-    "off_diag = np.sqrt(3.0/2.0) * torque_constant * gear_ratio\n",
-    "pto_impedance_12 = -1*(off_diag+0j) * np.ones(omega.shape) \n",
-    "pto_impedance_21 = -1*(off_diag+0j) * np.ones(omega.shape)\n",
-    "pto_impedance_22 = winding_impedance\n",
-    "pto_impedance_2 = np.array([[pto_impedance_11, pto_impedance_12],\n",
-    "                            [pto_impedance_21, pto_impedance_22]])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Next, we will create a new `PTO` object with this impedance matrix. We will also update the definitions of our PTO constraint and additional dynamic forcing function to use the new object. We will set our PTO constraint to $600 N$ for this example, since the dynamics for optimal electrical power will be different."
+    "# Optimal\n",
+    "x_wec = np.array([\n",
+    "         -5.34590357e-05,\n",
+    "        # -1.44175178e-04, -1.69577759e-04,\n",
+    "        # -1.93732274e-05,  2.95446742e-04,\n",
+    "         9.23568964e-04,  -3.34575436e-04,\n",
+    "        -6.84667669e-05, -4.50653781e-04,\n",
+    "         9.28453882e-04,   7.75556889e-05,\n",
+    "         1.65461770e-01,  -2.25669377e-01,\n",
+    "         5.34106551e-04,   1.36036042e-04,\n",
+    "        -1.36430591e-04,  4.02670813e-04,\n",
+    "        -3.80444249e-05,  1.63276367e-05,\n",
+    "        -1.57437036e-04,\n",
+    "        ]);\n",
+    "\n",
+    "x_opt = np.array([\n",
+    "        -1.30402770e+00,\n",
+    "        # -3.48316585e+00, -4.09789281e+00,\n",
+    "        # -4.67649366e-01,  6.92938769e+00,\n",
+    "         2.06482440e+01, -7.26660182e+00,\n",
+    "        -1.16610808e+00, -9.37008614e+00,\n",
+    "         1.74020260e+01,  2.51384525e+00,\n",
+    "         2.35013292e+03, -4.10306916e+03,\n",
+    "         7.42264452e+00,  3.45914541e+00,\n",
+    "        -3.14415012e+00,  4.46682294e+00,\n",
+    "        -4.50072846e-01, -6.55878884e-03,\n",
+    "        -1.16075834e+00,\n",
+    "        ]);"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [],
    "source": [
-    "## Update PTO\n",
-    "name_2 = ['PTO_Heave_Ex2']\n",
-    "pto_2 = wot.pto.PTO(ndof, kinematics, controller, pto_impedance_2, loss, name_2)\n",
+    "x = np.concatenate([x_wec, x_opt])\n",
     "\n",
-    "## Update PTO constraints and forcing\n",
-    "f_max_2 = 600.0\n",
-    "def const_f_pto_2(wec, x_wec, x_opt, waves):\n",
-    "    f = pto_2.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)\n",
-    "    return f_max_2 - np.abs(f.flatten())\n",
-    "ineq_cons_2 = {'type': 'ineq', 'fun': const_f_pto_2}\n",
-    "constraints_2 = [ineq_cons_2]\n",
-    "f_add_2 = {'PTO': pto_2.force_on_wec}"
+    "t = wec.time"
    ]
   },
   {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Finally, let's update our `WEC` object with the new PTO constraint, then run our optimization problem. Note we're now using `average_power` instead of `mechanical_average_power` as our objective function."
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "ValueError",
+     "evalue": "operands could not be broadcast together with shapes (32,) (36,) ",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[21], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# objetive function\u001b[39;00m\n\u001b[1;32m      2\u001b[0m f \u001b[38;5;241m=\u001b[39m problem[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfun\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[0;32m----> 3\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m      4\u001b[0m y\n",
+      "File \u001b[0;32m~/repos/WecOptTool/wecopttool/core.py:816\u001b[0m, in \u001b[0;36mWEC.solve.<locals>.obj_fun_scaled\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    815\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mobj_fun_scaled\u001b[39m(x):\n\u001b[0;32m--> 816\u001b[0m     x_wec, x_opt \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdecompose_state(\u001b[43mx\u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43mscale\u001b[49m)\n\u001b[1;32m    817\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m obj_fun(\u001b[38;5;28mself\u001b[39m, x_wec, x_opt, waves)\u001b[38;5;241m*\u001b[39mscale_obj\u001b[38;5;241m*\u001b[39msign\n",
+      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (32,) (36,) "
+     ]
+    }
+   ],
+   "source": [
+    "# objetive function\n",
+    "f = problem['fun']\n",
+    "y = f(x)\n",
+    "y"
    ]
   },
   {
@@ -510,43 +568,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Update WEC\n",
-    "\n",
-    "wec_2 = wot.WEC.from_bem(bem_data,\n",
-    "                         constraints=constraints_2,\n",
-    "                         friction=None,\n",
-    "                         f_add=f_add_2\n",
-    ")\n",
-    "\n",
-    "# Update objective function\n",
-    "obj_fun_2 = pto_2.average_power\n",
-    "\n",
-    "# Solve\n",
-    "scale_x_wec = 1e1  \n",
-    "scale_x_opt = 1e-3  \n",
-    "scale_obj = 1e-2  \n",
-    "\n",
-    "results_2 = wec_2.solve(\n",
-    "    waves, \n",
-    "    obj_fun_2, \n",
-    "    nstate_opt, \n",
-    "    scale_x_wec=scale_x_wec,\n",
-    "    scale_x_opt=scale_x_opt,\n",
-    "    scale_obj=scale_obj,\n",
-    ")\n",
-    "opt_average_power = results_2.fun\n",
-    "print(f'Optimal average electrical power: {opt_average_power} W')\n",
-    "\n",
-    "# Post-process\n",
-    "wec_fdom_2, wec_tdom_2 = wec_2.post_process(results_2, waves, nsubsteps)\n",
-    "pto_fdom_2, pto_tdom_2 = pto_2.post_process(wec_2, results_2, waves, nsubsteps)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We will compare our optimal results to the unconstrained case to gain some insight into the effect of the constraint on the optimal PTO force. Let's do the same process as before, but unset the `constraints` parameter in a new `WEC` object."
+    "# jacobian\n",
+    "f = problem['jac']\n",
+    "y = f(x)\n",
+    "plt.plot(y, 'o')"
    ]
   },
   {
@@ -555,34 +580,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "wec_2_nocon = wot.WEC.from_bem(\n",
-    "    bem_data,\n",
-    "    constraints=None,\n",
-    "    friction=None,\n",
-    "    f_add=f_add_2)\n",
-    "\n",
-    "results_2_nocon = wec_2_nocon.solve(\n",
-    "    waves,\n",
-    "    obj_fun_2, \n",
-    "    nstate_opt,\n",
-    "    scale_x_wec=scale_x_wec,\n",
-    "    scale_x_opt=scale_x_opt,\n",
-    "    scale_obj=scale_obj,\n",
-    ")\n",
-    "opt_average_power = results_2_nocon.fun\n",
-    "print(f'Optimal average electrical power: {opt_average_power} W')\n",
-    "wec_fdom_2_nocon, wec_tdom_2_nocon = wec_2_nocon.post_process(\n",
-    "    results_2_nocon, waves, nsubsteps)\n",
-    "pto_fdom_2_nocon, pto_tdom_2_nocon = pto_2.post_process(\n",
-    "    wec_2_nocon, results_2_nocon, waves, nsubsteps)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Note that the optimal constrained PTO force follows the optimal unconstrained solution (sinusoidal) whenever the unconstrained solution is within the constraint. \n",
-    "When the constraint is active the optimal PTO force is the maximum PTO force of $600 N$."
+    "# constraint (dynamics in residual form)\n",
+    "problem['constraints']"
    ]
   },
   {
@@ -591,59 +590,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.figure()\n",
-    "wec_tdom_2['pos'].plot(label='constrained')\n",
-    "wec_tdom_2_nocon['pos'].plot(label='unconstrained')\n",
-    "plt.legend(loc='lower right')\n",
-    "\n",
-    "plt.figure()\n",
-    "pto_tdom_2['force'].plot(label='constrained')\n",
-    "pto_tdom_2_nocon['force'].plot(label='unconstrained')\n",
-    "plt.legend(loc='lower right')\n",
-    "\n",
-    "plt.figure()\n",
-    "pto_tdom_2['power'].loc['elec',:,:].plot(label='constrained')\n",
-    "pto_tdom_2_nocon['power'].loc['elec',:,:].plot(label='unconstrained')\n",
-    "plt.legend(loc='lower right')\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The attentive user might have noticed that the amplitude of the position, force and power signals is about half the magnitude of the signals we plotted in the first part of the tutorial. We can see that optimizing for electrical power requires optimal state trajectories with smaller amplitudes. For most WECs the electrical power is the usable form of power, thus the WEC should be designed for electrical power and we can avoid over-designing, which would results from expecting the forces associated with the optimal trajectories for mechanical power maximisation."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## 3. Control co-design of the WEC geometry for maximum electrical power\n",
-    "The first two examples only used the inner optimization loop in WecOptTool to optimize PTO power. Here in Part 3 we bring it all together and show how to use both the inner and outer optimization loops in WecOptTool to do control co-optimization of a hull design in conjunction with an optimal controller for electrical power.\n",
-    "Again, we use the WaveBot WEC in one degree of freedom in regular waves. \n",
-    "The goal is to **find the optimal keel radius** (`r2`) that maximizes the average produced electrical power, while maintaining a constant hull volume. \n",
-    "A constant volume is achieved by setting the height of the conical section (`h2`) in conjunction with the keel radius (`r2`).\n",
-    "\n",
-    "This example demonstrates a complete case of the types of optimization studies WecOptTool is meant for. \n",
-    "The main optimization (outer optimization loop) is to find the optimal geometry (radius `r2`), and for each geometry considered the optimal PTO force (inner optimization loop) will be found.\n",
-    "The inner loop was showcased in Example 2 and uses a gradient-based optimization method, with the gradients obtained with automatic differentiation. \n",
-    "The outer loop optimization is for the user to setup. \n",
-    "In this example, we will do a simple *brute force* optimization using `scipy.optimize.brute`.  \n",
-    "\n",
-    "![Device Diagram](https://live.staticflickr.com/65535/51751577441_515afec334_z.jpg) \n",
-    "<div>\n",
-    "<img src=\"https://live.staticflickr.com/65535/52434071157_187eb4334c_k.jpg\" width=\"1000\">\n",
-    "</div>\n",
-    "\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Problem setup\n",
-    "First, we define a function for `h2` based on `r1` that maintains a constant volume. \n",
-    "We see that, as expected, smaller values of `r2` require larger values of `h2` in order to maintain a constant hull volume."
+    "f = problem['constraints'][0]['fun']\n",
+    "y = f(x)\n",
+    "plt.plot(t, y)"
    ]
   },
   {
@@ -652,36 +601,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "r1 = 0.88\n",
-    "r2_0 = 0.35\n",
-    "h2_0 = 0.37\n",
-    "V0 = 1/3*np.pi*h2_0*(r1**2+r2_0**2+(r1*r2_0))\n",
-    "\n",
-    "r2_vals = np.linspace(0.05, 0.88*0.999, 8, endpoint=True)\n",
-    "\n",
-    "\n",
-    "def h2_from_r2(r2, V=V0, r1=r1):\n",
-    "    h2 = V/(1/3*np.pi*(r1**2+r2**2+(r1*r2)))\n",
-    "    return h2\n",
-    "\n",
-    "\n",
-    "# plot\n",
-    "mapres = map(h2_from_r2, r2_vals)\n",
-    "h2_vals = list(mapres)\n",
-    "\n",
-    "fig1, ax1 = plt.subplots(figsize=(8,5))\n",
-    "for r2, h2 in zip(r2_vals.tolist(), h2_vals):\n",
-    "    _ = wot.geom.WaveBot(r2=r2, h2=h2, freeboard=0.2).plot_cross_section(\n",
-    "        ax=ax1, label=f\"r2={r2:.2f}, h2={h2:.2f}\")\n",
-    "ax1.legend(loc='best', fontsize='small',ncol=2)\n",
-    "_ = ax1.set_title('WaveBot hull cross-sections')\n"
+    "f = problem['constraints'][0]['jac']\n",
+    "y = f(x)\n",
+    "plt.imshow(y, cmap=\"seismic\")\n",
+    "plt.colorbar()\n",
+    "print(y.shape)"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "Next we will define an objective function for our design optimization problem. We use the same workflow illustrated in Part 2 to set up a WaveBot device and solve for the optimal solution, but wrap this in a function definition which can set `r2` and (indirectly) `h2`."
+    "f = problem['constraints'][0]['jac']\n",
+    "y = f(x) * x\n",
+    "plt.imshow(y, cmap=\"seismic\")\n",
+    "plt.colorbar()"
    ]
   },
   {
@@ -690,114 +626,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def design_obj_fun(x):\n",
-    "\n",
-    "    # Unpack geometry variables\n",
-    "    r2 = x[0]\n",
-    "    h2 = h2_from_r2(r2)\n",
-    "    print(f\"\\nr2 = {r2:.2f}:\")\n",
-    "\n",
-    "    # Set up Capytaine floating body\n",
-    "    wb = wot.geom.WaveBot(r2=r2, h2=h2)\n",
-    "    mesh = wb.mesh(mesh_size_factor=0.5)\n",
-    "    fb = cpy.FloatingBody.from_meshio(mesh, name=\"WaveBot\")\n",
-    "    fb.add_translation_dof(name=\"Heave\")\n",
-    "    ndof = fb.nb_dofs\n",
-    "\n",
-    "    # Run BEM\n",
-    "    f1 = 0.05\n",
-    "    nfreq = 50\n",
-    "    bem_data = wot.run_bem(fb, freq)\n",
-    "\n",
-    "    # Impedance definition\n",
-    "    omega = bem_data.omega.values\n",
-    "    gear_ratio = 12.0\n",
-    "    torque_constant = 6.7\n",
-    "    winding_resistance = 0.5\n",
-    "    winding_inductance = 0.0\n",
-    "    drivetrain_inertia = 2.0\n",
-    "    drivetrain_friction = 1.0\n",
-    "    drivetrain_stiffness = 0.0\n",
-    "\n",
-    "    drivetrain_impedance = (1j*omega*drivetrain_inertia + \n",
-    "                            drivetrain_friction + \n",
-    "                            1/(1j*omega)*drivetrain_stiffness) \n",
-    "\n",
-    "    winding_impedance = winding_resistance + 1j*omega*winding_inductance\n",
-    "\n",
-    "\n",
-    "    pto_impedance_11 = -1* gear_ratio**2 * drivetrain_impedance\n",
-    "    off_diag = np.sqrt(3.0/2.0) * torque_constant * gear_ratio\n",
-    "    pto_impedance_12 = -1*(off_diag+0j) * np.ones(omega.shape) \n",
-    "    pto_impedance_21 = -1*(off_diag+0j) * np.ones(omega.shape)\n",
-    "    pto_impedance_22 = winding_impedance\n",
-    "    pto_impedance_3 = np.array([[pto_impedance_11, pto_impedance_12],\n",
-    "                            [pto_impedance_21, pto_impedance_22]])\n",
-    "\n",
-    "    # Set PTO object\n",
-    "    name = [\"PTO_Heave\",]\n",
-    "    kinematics = np.eye(ndof)\n",
-    "    efficiency = None\n",
-    "    controller = None\n",
-    "    pto = wot.pto.PTO(ndof, kinematics, controller, pto_impedance_3, efficiency, name)\n",
-    "    \n",
-    "\n",
-    "    # Set PTO constraint and additional dynamic force\n",
-    "    nsubsteps = 4\n",
-    "    f_max = 600.0\n",
-    "\n",
-    "    def const_f_pto(wec, x_wec, x_opt, waves):\n",
-    "        f = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)\n",
-    "        return f_max - np.abs(f.flatten())\n",
-    "\n",
-    "    ineq_cons = {'type': 'ineq', 'fun': const_f_pto}\n",
-    "    constraints = [ineq_cons]\n",
-    "\n",
-    "    f_add = {'PTO': pto.force_on_wec}\n",
-    "\n",
-    "    # Create WEC\n",
-    "    wec = wot.WEC.from_bem(bem_data,\n",
-    "            constraints=constraints,\n",
-    "            friction=None, \n",
-    "            f_add=f_add,\n",
-    "            )\n",
-    "\n",
-    "    # Waves\n",
-    "    wfreq = 0.3\n",
-    "    amplitude = 0.0625\n",
-    "    phase = -40\n",
-    "    waves_3 = wot.waves.regular_wave(f1, nfreq, wfreq, amplitude, phase)\n",
-    "\n",
-    "    # Objective function\n",
-    "    obj_fun = pto.average_power\n",
-    "    nstate_opt = 2*nfreq\n",
-    "    \n",
-    "    # Solve\n",
-    "    scale_x_wec = 1e1  \n",
-    "    scale_x_opt = 1e-3  \n",
-    "    scale_obj = 1e-2  \n",
-    "    res = wec.solve(\n",
-    "        waves, \n",
-    "        obj_fun, \n",
-    "        nstate_opt, \n",
-    "        scale_x_wec=scale_x_wec,\n",
-    "        scale_x_opt=scale_x_opt,\n",
-    "        scale_obj=scale_obj)\n",
-    "\n",
-    "    return res.fun"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Solve\n",
-    "Finally, we may call this objective function with an optimization algorithm. \n",
-    "Here, a simple *brute force* optimization approach is used for illustrative purposes, but any variety of options could be applied. \n",
-    "The optimization algorithm will call our objective function, which in turn will create a new WaveBot hull, run the necessary BEM calculations for the hull, and find the PTO force that provides the most electric power for that hull. \n",
-    "This process will be conducted for the range of `r2` values that we specify.\n",
-    "\n",
-    "_(note: the cell below will likely take 5+ minutes to run on a standard personal computer)_"
+    "# callback\n",
+    "f = problem['callback']\n",
+    "y = f(x)\n",
+    "y == None"
    ]
   },
   {
@@ -806,47 +638,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "wot.set_loglevel(\"error\")  # Suppress warnings\n",
-    "\n",
-    "# range over which to search\n",
-    "ranges = (slice(r2_vals[0], r2_vals[-1]+np.diff(r2_vals)[0], np.diff(r2_vals)[0]),)\n",
-    "\n",
-    "# solve\n",
-    "opt_x0, opt_fval, x0s, fvals = brute(func=design_obj_fun, ranges=ranges, full_output=True,  finish=None)"
+    "# x0 initial guess\n",
+    "plt.plot(problem['x0'], 'o')\n",
+    "plt.plot([16, 16], [-3, 3], 'k--')"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
-   "source": [
-    "### Results\n",
-    "From a quick plot of the results, we see that the power absorption (where negative power is power absorbed by the device) generally improves for smaller values of `r2`.\n",
-    "It is also clear that when the WEC is cylindrical (where `r2=0.88`), power absorption is reduced."
-   ]
+   "outputs": [],
+   "source": []
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": [
-    "fig, ax = plt.subplots()\n",
-    "colors = plt.rcParams['axes.prop_cycle'].by_key()['color'][:len(x0s)]\n",
-    "ax.plot(x0s, fvals, 'k', zorder=0)\n",
-    "ax.scatter(x0s, fvals, c=colors, zorder=1)\n",
-    "\n",
-    "ax.set_xlabel('Keel radius, $r_2$ [m]')\n",
-    "ax.set_ylabel('Average Power [W]')\n",
-    "ax.set_title('Design optimization results')\n",
-    "fig.tight_layout()"
-   ]
+   "source": []
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
-   "source": [
-    "Note that in this case the magnitude of average power between the different keel radii is rather small, this is because the PTO force constraint is active most of the time, therefore all considered geometries perform similarily. If you remove the PTO constraint and re-run the co-optimization study you will see that the impact of radius on average electrical power is significantly higher."
-   ]
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -865,7 +681,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.10.13"
   },
   "vscode": {
    "interpreter": {
diff --git a/wecopttool/core.py b/wecopttool/core.py
index a569c0da..6077382a 100644
--- a/wecopttool/core.py
+++ b/wecopttool/core.py
@@ -36,7 +36,7 @@
     "force_from_rao_transfer_function",
     "force_from_impedance",
     "force_from_waves",
-    "inertia",
+    # "inertia",
     "standard_forces",
     "run_bem",
     "change_bem_convention",
@@ -51,7 +51,7 @@
     "frequency_parameters",
     "time_results",
     "set_fb_centers",
-]
+    ]
 
 
 import logging
@@ -114,9 +114,10 @@ class WEC:
     """
     def __init__(
         self,
-        f1: float,
-        nfreq: int,
+        df: float,
+        ifreq_end: int,
         forces: TIForceDict,
+        ifreq_start: int = 1,
         constraints: Optional[Iterable[Mapping]] = None,
         inertia_matrix: Optional[ndarray] = None,
         ndof: Optional[int] = None,
@@ -206,11 +207,15 @@ def __init__(
             Initialize a :py:class:`wecopttool.WEC` object from an
             intrinsic impedance array and excitation coefficients.
         """
-        self._freq = frequency(f1, nfreq)
-        self._time = time(f1, nfreq)
-        self._time_mat = time_mat(f1, nfreq)
-        self._derivative_mat = derivative_mat(f1, nfreq)
-        self._derivative2_mat = derivative2_mat(f1, nfreq)
+        self._df = df
+        self._ifreq_start = ifreq_start
+        self._ifreq_end = ifreq_end
+        idx_1, idx_2 = self._idxs(ifreq_start, ifreq_end)
+        self._freq = frequency(df, ifreq_end)[idx_1]
+        self._time = time(df, ifreq_end)
+        self._time_mat = time_mat(df, ifreq_end)[:, idx_2]
+        self._derivative_mat = derivative_mat(df, ifreq_end)[idx_2][:, idx_2]
+        self._derivative2_mat = derivative2_mat(df, ifreq_end)[idx_2][:, idx_2]
         self._forces = forces
         constraints = list(constraints) if (constraints is not None) else []
         self._constraints = constraints
@@ -258,7 +263,7 @@ def _ignored(var_name, condition):
         if inertia_in_forces:
             _inertia = None
         else:
-            _inertia = inertia(f1, nfreq, inertia_matrix)
+            _inertia = WEC._get_inertia(df, ifreq_start, ifreq_end, inertia_matrix)
         self._inertia = _inertia
 
         # names
@@ -271,7 +276,7 @@ def _ignored(var_name, condition):
     def __str__(self) -> str:
         str = (f'{self.__class__.__name__}: ' +
                f'DOFs ({self.ndof})={self.dof_names}, ' +
-               f'f=[0, {self.f1}, ..., {self.nfreq}({self.f1})] Hz.')
+               f'f=[0, {self.df}, ..., {self.ifreq_end}({self.df})] Hz.')
         return str
 
     def __repr__(self) -> str:
@@ -281,6 +286,84 @@ def __repr__(self) -> str:
         repr_org = f"<{module}.{qualname} object at {hex(id(self))}>"
         return repr_org + " :: " + self.__str__()
 
+    @staticmethod
+    def _idxs(ifreq_start, ifreq_end):
+        idx_1 = [True if i==0 or i>=ifreq_start else False for i in range(ifreq_end+1)]
+        idx_2 = [True if i==0 or i>=2*ifreq_start-1 else False for i in range(2*ifreq_end)]
+        return idx_1, idx_2
+
+    @staticmethod
+    def _get_inertia(
+        df: float,
+        ifreq_start: int,
+        ifreq_end: int,
+        inertia_matrix: ArrayLike,
+        ) -> TStateFunction:
+        """Create the inertia "force" from the inertia matrix.
+
+        Parameters
+        ----------
+        df
+            Frequency spacing [:math:`Hz`].
+        nfreq
+            Number of frequencies.
+        inertia_matrix
+            Inertia matrix.
+        """
+        omega = np.expand_dims(frequency(df, ifreq_end)[ifreq_start:]*2*np.pi, [1,2])
+        inertia_matrix = np.expand_dims(inertia_matrix, 0)
+        rao_transfer_function = -1*omega**2*inertia_matrix + 0j
+        inertia_fun = force_from_rao_transfer_function(
+            rao_transfer_function, False)
+        return inertia_fun
+
+    @staticmethod
+    def frequency_parameters(
+        freqs: ArrayLike,
+        precision: Optional[int] = 10,
+        ) -> tuple[float, int]:
+        """Return the fundamental frequency, the number of frequencies,
+        and first frequency in a frequency array.
+
+        This function can be used as a check for inputs to other
+        functions since it raises an error if the frequency vector does
+        not have the correct format (equally spaced).
+
+        Parameters
+        ----------
+        freqs
+            The frequency array with equal spacing.
+        precision
+            Controls rounding of fundamental frequency.
+
+        Returns
+        -------
+        df
+            Frequency spacing [:math:`Hz`].
+        ifreq_end
+            Last frequency (index).
+        ifreq_start
+            Frequency (index) at which vector starts.
+
+        Raises
+        ------
+        ValueError
+            If the frequency vector is not evenly spaced.
+        """
+        df = freqs[1] - freqs[0]
+        df = df if precision is None else round(df, precision)
+        ifreq_start = round(freqs[0]/df)
+        assert np.isclose(ifreq_start, freqs[0]/df)
+        ifreq_end = ifreq_start + len(freqs) - 1
+        f_check = np.arange(ifreq_start, ifreq_end+1)*df
+        if not np.allclose(f_check, freqs):
+            print(f_check)
+            print(freqs)
+            raise ValueError(
+                "Frequency array must be evenly spaced by " +
+                "the fundamental frequency ")
+        return df, ifreq_start, ifreq_end
+
     # other initialization methods
     @staticmethod
     def from_bem(
@@ -311,8 +394,8 @@ def from_bem(
         :py:func:`wecopttool.run_bem`,
         rather than running Capytaine directly, which outputs the
         results in the correct convention. The results can be saved
-        using :py:func:`wecopttool.write_netcdf`. 
-        :py:func:`wecopttool.run_bem` also computes the inertia and 
+        using :py:func:`wecopttool.write_netcdf`.
+        :py:func:`wecopttool.run_bem` also computes the inertia and
         hydrostatic stiffness which should be included in bem_data.
 
         Parameters
@@ -363,8 +446,8 @@ def from_bem(
         inertia_matrix = hydro_data['inertia_matrix'].values
 
         # frequency array
-        f1, nfreq = frequency_parameters(
-            hydro_data.omega.values/(2*np.pi), False)
+        df, ifreq_start, ifreq_end = WEC.frequency_parameters(
+            hydro_data.omega.values/(2*np.pi))
 
         # check real part of damping diagonal > 0
         if min_damping is not None:
@@ -377,13 +460,17 @@ def from_bem(
         forces = linear_force_functions | f_add
         # constraints
         constraints = constraints if (constraints is not None) else []
-        return WEC(f1, nfreq, forces, constraints, inertia_matrix, dof_names=dof_names)
+        wec = WEC(
+            df, ifreq_end, forces, ifreq_start, constraints,
+            inertia_matrix, dof_names=dof_names)
+        return wec
 
     @staticmethod
     def from_floating_body(
         fb: cpy.FloatingBody,
-        f1: float,
-        nfreq: int,
+        df: float,
+        ifreq_end: int,
+        ifreq_start: int = 1,
         friction: Optional[ndarray] = None,
         f_add: Optional[TIForceDict] = None,
         constraints: Optional[Iterable[Mapping]] = None,
@@ -392,7 +479,7 @@ def from_floating_body(
         rho: Optional[float] = _default_parameters['rho'],
         g: Optional[float] = _default_parameters['g'],
         depth: Optional[float] = _default_parameters['depth'],
-    ) -> TWEC:
+        ) -> TWEC:
         """Create a WEC object from a Capytaine :python:`FloatingBody`
         (:py:class:capytaine.bodies.bodies.FloatingBody).
 
@@ -464,9 +551,11 @@ def from_floating_body(
         """
 
         # RUN BEM
-        _log.info(f"Running Capytaine (BEM): {nfreq+1} frequencies x " +
+        nfreq = ifreq_end - ifreq_start + 1
+        _log.info(f"Running Capytaine (BEM): {nfreq} frequencies x " +
                  f"{len(wave_directions)} wave directions.")
-        freq = frequency(f1, nfreq)[1:]
+        idx_1, _ = self._idxs(ifreq_start, ifreq_end)
+        freq = frequency(df, ifreq_end)[idx_1][1:]
         bem_data = run_bem(
             fb, freq, wave_directions, rho=rho, g=g, depth=depth)
         wec = WEC.from_bem(
@@ -483,7 +572,7 @@ def from_impedance(
         f_add: Optional[TIForceDict] = None,
         constraints: Optional[Iterable[Mapping]] = None,
         min_damping: Optional[float] = _default_min_damping,
-    ) -> TWEC:
+        ) -> TWEC:
         """Create a WEC object from the intrinsic impedance and
         excitation coefficients.
 
@@ -496,7 +585,6 @@ def from_impedance(
         The impedance can also be obtained experimentally.
         Note that the impedance is not defined at :math:`ω=0`.
 
-
         Parameters
         ----------
         freqs
@@ -532,7 +620,8 @@ def from_impedance(
             If :python:`impedance` does not have the correct size:
             :python:`(ndof, ndof, nfreq)`.
         """
-        f1, nfreq = frequency_parameters(freqs, False)
+        df, ifreq_start, ifreq_end = self.frequency_parameters(freqs)
+        nfreq = ifreq_end - ifreq_start + 1
 
         # impedance matrix shape
         shape = impedance.shape
@@ -564,10 +653,11 @@ def from_impedance(
         forces = forces | f_add
 
         # wec
-        wec = WEC(f1, nfreq, forces, constraints,
-                  inertia_in_forces=True, ndof=shape[1])
+        wec = WEC(df, ifreq_end, forces, ifreq_start, constraints,
+            inertia_in_forces=True, ndof=shape[1])
         return wec
 
+    # solve
     def residual(self, x_wec: ndarray, x_opt: ndarray, waves: Dataset,
         ) -> float:
         """
@@ -586,13 +676,12 @@ def residual(self, x_wec: ndarray, x_opt: ndarray, waves: Dataset,
         if not self.inertia_in_forces:
             ri = self.inertia(self, x_wec, x_opt, waves)
         else:
-            ri = np.zeros([self.ncomponents, self.ndof])
+            ri = np.zeros([self.nt, self.ndof])
         # forces, -Σf
         for f in self.forces.values():
             ri = ri - f(self, x_wec, x_opt, waves)
         return self.dofmat_to_vec(ri)
 
-    # solve
     def solve(self,
         waves: Dataset,
         obj_fun: TStateFunction,
@@ -795,17 +884,19 @@ def callback_scipy(x):
 
         # optimization problem
         optim_options['disp'] = optim_options.get('disp', True)
-        problem = {'fun': obj_fun_scaled,
-                    'x0': x0,
-                    'method': 'SLSQP',
-                    'constraints': constraints,
-                    'options': optim_options,
-                    'bounds': bounds,
-                    'callback': callback_scipy,
-                    }
+        problem = {
+            'fun': obj_fun_scaled,
+            'x0': x0,
+            'method': 'SLSQP',
+            'constraints': constraints,
+            'options': optim_options,
+            'bounds': bounds,
+            'callback': callback_scipy,
+        }
         if use_grad:
             problem['jac'] = grad(obj_fun_scaled)
 
+        return problem ## DEBUG
         # minimize
         optim_res = minimize(**problem)
 
@@ -828,7 +919,7 @@ def post_process(self,
         res: OptimizeResult,
         waves: Dataset,
         nsubsteps: Optional[int] = 1,
-    ) -> tuple[Dataset, Dataset]:
+        ) -> tuple[Dataset, Dataset]:
         """Post-process the results from
         :py:meth:`wecopttool.WEC.solve`.
 
@@ -1006,15 +1097,30 @@ def frequency(self) -> ndarray:
         """Frequency vector [:math:`Hz`]."""
         return self._freq
 
+    @property
+    def df(self) -> float:
+        """Frequency spacing [:math:`Hz`]."""
+        return self._df
+
     @property
     def f1(self) -> float:
         """Fundamental frequency :python:`f1` [:math:`Hz`]."""
-        return self._freq[1]
+        return self.df
+
+    @property
+    def ifreq_end(self) -> int:
+        """Number of frequencies, not including the zero-frequency."""
+        return self._ifreq_end
+
+    @property
+    def ifreq_start(self) -> int:
+        """Number of frequencies, not including the zero-frequency."""
+        return self._ifreq_start
 
     @property
     def nfreq(self) -> int:
         """Number of frequencies, not including the zero-frequency."""
-        return len(self._freq)-1
+        return self.ifreq_end - self.ifreq_start + 1
 
     @property
     def omega(self) -> ndarray:
@@ -1085,12 +1191,12 @@ def tf(self) -> float:
         """Final time (repeat period) [s]. Not included in
         :python:`time` vector.
         """
-        return 1/self.f1
+        return 1/self.df
 
     @property
     def nt(self) -> int:
         """Number of timesteps."""
-        return self.ncomponents
+        return 2*self.ifreq_end
 
     @property
     def ncomponents(self) -> int:
@@ -1111,7 +1217,7 @@ def nstate_wec(self) -> int:
     # other methods
     def decompose_state(self,
         state: ndarray
-    ) -> tuple[ndarray, ndarray]:
+        ) -> tuple[ndarray, ndarray]:
         """Split the state vector into the WEC dynamics state and the
         optimization (control) state.
 
@@ -1155,7 +1261,7 @@ def time_nsubsteps(self, nsubsteps: int) -> ndarray:
         --------
         time, WEC.time
         """
-        return time(self.f1, self.nfreq, nsubsteps)
+        return time(self.df, self.ifreq_end, nsubsteps)
 
     def time_mat_nsubsteps(self, nsubsteps: int) -> ndarray:
         """Create a time matrix similar to
@@ -1174,7 +1280,8 @@ def time_mat_nsubsteps(self, nsubsteps: int) -> ndarray:
         --------
         time_mat, WEC.time_mat, WEC.time_nsubsteps
         """
-        return time_mat(self.f1, self.nfreq, nsubsteps)
+        _, idx_2 = self._idxs(self.ifreq_start, self.ifreq_end)
+        return time_mat(self.df, self.ifreq_end, nsubsteps)[:, idx_2]
 
     def vec_to_dofmat(self, vec: ndarray) -> ndarray:
         """Convert a vector to a matrix with one column per degree of
@@ -1238,7 +1345,7 @@ def fd_to_td(self, fd: ndarray) -> ndarray:
         --------
         fd_to_td, WEC.td_to_fd
         """
-        return fd_to_td(fd, self.f1, self.nfreq, True)
+        return np.dot(self.time_mat, complex_to_real(fd))
 
     def td_to_fd(
         self,
@@ -1264,13 +1371,14 @@ def td_to_fd(
         --------
         td_to_fd, WEC.fd_to_td
         """
-        return td_to_fd(td, fft, True)
+        idx_1, _ = self._idxs(self.ifreq_start, self.ifreq_end)
+        return td_to_fd(td, fft, True)[idx_1, :]
 
 
 def ncomponents(
     nfreq : int,
     zero_freq: Optional[bool] = True,
-) -> int:
+    ) -> int:
     """Number of Fourier components (:python:`2*nfreq`) for each
     DOF. The sine component of the highest frequency (the 2-point wave)
     is excluded as it will always evaluate to zero.
@@ -1295,7 +1403,7 @@ def frequency(
     f1: float,
     nfreq: int,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Construct equally spaced frequency array.
 
     The array includes :python:`0` and has length of :python:`nfreq+1`.
@@ -1325,7 +1433,7 @@ def time(
     f1: float,
     nfreq: int,
     nsubsteps: Optional[int] = 1,
-) -> ndarray:
+    ) -> ndarray:
     """Assemble the time vector with :python:`nsubsteps` subdivisions.
 
     Returns the 1D time vector, in seconds, starting at time
@@ -1355,7 +1463,7 @@ def time_mat(
     nfreq: int,
     nsubsteps: Optional[int] = 1,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Assemble the time matrix that converts the state to a
     time-series.
 
@@ -1400,7 +1508,7 @@ def derivative_mat(
     f1: float,
     nfreq: int,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Assemble the derivative matrix that converts the state vector of
     a response to the state vector of its derivative.
 
@@ -1437,7 +1545,7 @@ def derivative2_mat(
     f1: float,
     nfreq: int,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Assemble the second derivative matrix that converts the state vector of
     a response to the state vector of its second derivative.
 
@@ -1472,7 +1580,7 @@ def derivative2_mat(
 def mimo_transfer_mat(
     transfer_mat: DataArray,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Create a block matrix of the MIMO transfer function.
 
     The input is a complex transfer matrix that relates the complex
@@ -1568,7 +1676,7 @@ def dofmat_to_vec(mat: ArrayLike) -> ndarray:
 def real_to_complex(
     fd: ArrayLike,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Convert from two real amplitudes to one complex amplitude per
     frequency.
 
@@ -1617,7 +1725,7 @@ def real_to_complex(
 def complex_to_real(
     fd: ArrayLike,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Convert from one complex amplitude to two real amplitudes per
     frequency.
 
@@ -1679,7 +1787,7 @@ def fd_to_td(
     f1: Optional[float] = None,
     nfreq: Optional[int] = None,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Convert a complex array of Fourier coefficients to a real array
     of time-domain responses.
 
@@ -1748,7 +1856,7 @@ def td_to_fd(
     td: ArrayLike,
     fft: Optional[bool] = True,
     zero_freq: Optional[bool] = True,
-) -> ndarray:
+    ) -> ndarray:
     """Convert a real array of time-domain responses to a complex array
     of Fourier coefficients.
 
@@ -1828,7 +1936,7 @@ def check_linear_damping(
     hydro_data: Dataset,
     min_damping: Optional[float] = 1e-6,
     uniform_shift: Optional[bool] = True,
-) -> Dataset:
+    ) -> Dataset:
     """Ensure that the linear hydrodynamics (friction + radiation
     damping) have positive damping.
 
@@ -1885,7 +1993,7 @@ def check_linear_damping(
 def check_impedance(
     Zi: DataArray,
     min_damping: Optional[float] = 1e-6,
-) -> DataArray:
+    ) -> DataArray:
     """Ensure that the real part of the impedance (resistive) is positive.
 
     Adds to real part of the impedance.
@@ -1916,7 +2024,7 @@ def check_impedance(
 def force_from_rao_transfer_function(
     rao_transfer_mat: DataArray,
     zero_freq: Optional[bool] = True,
-) -> TStateFunction:
+    ) -> TStateFunction:
     """Create a force function from its position transfer matrix.
 
     This is the position equivalent to the velocity-based
@@ -1947,7 +2055,7 @@ def force(wec, x_wec, x_opt, waves):
 def force_from_impedance(
     omega: ArrayLike,
     impedance: DataArray,
-) -> TStateFunction:
+    ) -> TStateFunction:
     """Create a force function from its impedance.
 
     Parameters
@@ -1980,30 +2088,6 @@ def force(wec, x_wec, x_opt, waves):
     return force
 
 
-def inertia(
-    f1: float,
-    nfreq: int,
-    inertia_matrix: ArrayLike,
-) -> TStateFunction:
-    """Create the inertia "force" from the inertia matrix.
-
-    Parameters
-    ----------
-    f1
-        Fundamental frequency :python:`f1` [:math:`Hz`].
-    nfreq
-        Number of frequencies.
-    inertia_matrix
-        Inertia matrix.
-    """
-    omega = np.expand_dims(frequency(f1, nfreq, False)*2*np.pi, [1,2])
-    inertia_matrix = np.expand_dims(inertia_matrix, 0)
-    rao_transfer_function = -1*omega**2*inertia_matrix + 0j
-    inertia_fun = force_from_rao_transfer_function(
-        rao_transfer_function, False)
-    return inertia_fun
-
-
 def standard_forces(hydro_data: Dataset) -> TForceDict:
     """Create functions for linear hydrodynamic forces.
 
@@ -2064,7 +2148,7 @@ def run_bem(
     depth: float = _default_parameters['depth'],
     write_info: Optional[Mapping[str, bool]] = None,
     njobs: int = 1,
-) -> Dataset:
+    ) -> Dataset:
     """Run Capytaine for a range of frequencies and wave directions.
 
     This simplifies running *Capytaine* and ensures the output are in
@@ -2160,7 +2244,7 @@ def change_bem_convention(bem_data: Dataset) -> Dataset:
 def add_linear_friction(
     bem_data: Dataset,
     friction: Optional[ArrayLike] = None
-) -> Dataset:
+    ) -> Dataset:
     """Add linear friction to BEM data.
 
     Returns the Dataset with the additional information added.
@@ -2178,7 +2262,7 @@ def add_linear_friction(
     """
     dims = ['radiating_dof', 'influenced_dof']
     hydro_data = bem_data.copy(deep=True)
-    
+
     if friction is not None:
         if 'friction' in hydro_data.variables.keys():
             if not np.allclose(data, hydro_data.variables[name]):
@@ -2195,7 +2279,7 @@ def add_linear_friction(
     elif friction is None:
         ndof = len(hydro_data["influenced_dof"])
         hydro_data['friction'] = (dims, np.zeros([ndof, ndof]))
-            
+
     return hydro_data
 
 
@@ -2286,7 +2370,7 @@ def atleast_2d(array: ArrayLike) -> ndarray:
 def degrees_to_radians(
     degrees: FloatOrArray,
     sort: Optional[bool] = True,
-) -> Union[float, ndarray]:
+    ) -> Union[float, ndarray]:
     """Convert a 1D array of angles in degrees to radians in the range
     :math:`[-π, π)` and optionally sort them.
 
@@ -2313,7 +2397,7 @@ def subset_close(
     rtol: float = 1.e-5,
     atol: float = 1.e-8,
     equal_nan: bool = False,
-) -> tuple[bool, list]:
+    ) -> tuple[bool, list]:
     """Check if the first set :python:`set_a` is contained, to some
     tolerance, in the second set :python:`set_b`.
 
@@ -2400,7 +2484,7 @@ def decompose_state(
     state: ndarray,
     ndof: int,
     nfreq: int,
-) -> tuple[ndarray, ndarray]:
+    ) -> tuple[ndarray, ndarray]:
     """Split the state vector into the WEC dynamics state and the
     optimization (control) state.
 
@@ -2432,7 +2516,7 @@ def decompose_state(
 def frequency_parameters(
     freqs: ArrayLike,
     zero_freq: bool = True,
-) -> tuple[float, int]:
+    ) -> tuple[float, int]:
     """Return the fundamental frequency and the number of frequencies
     in a frequency array.
 
@@ -2508,13 +2592,13 @@ def time_results(fd: DataArray, time: DataArray) -> ndarray:
 def set_fb_centers(
     fb: FloatingBody,
     rho: float = _default_parameters["rho"],
-) -> FloatingBody:
+    ) -> FloatingBody:
     """Sets default properties if not provided by the user:
         - `center_of_mass` is set to the geometric centroid
         - `rotation_center` is set to the center of mass
     """
     valid_properties = ['center_of_mass', 'rotation_center']
-    
+
     for property in valid_properties:
         if not hasattr(fb, property):
             setattr(fb, property, None)
@@ -2533,5 +2617,5 @@ def set_fb_centers(
         elif getattr(fb, property) is not None:
             _log.warning(
                 f'{property} already defined as {getattr(fb, property)}.')
-            
+
     return fb
diff --git a/wecopttool/pto.py b/wecopttool/pto.py
index e34b3b1e..a178560b 100644
--- a/wecopttool/pto.py
+++ b/wecopttool/pto.py
@@ -1080,7 +1080,11 @@ def controller_p(
 
 
 # utilities
-def nstate_unstructured(nfreq: int, ndof: int) -> int:
+def nstate_unstructured(
+        ndof: int,
+        ifreq_end: int,
+        ifreq_start=1
+    ) -> int:
     """
     Number of states needed to represent an unstructured controller.
 
@@ -1091,6 +1095,7 @@ def nstate_unstructured(nfreq: int, ndof: int) -> int:
     ndof
         Number of degrees of freedom.
     """
+    nfreq = ifreq_end - ifreq_start + 1
     return 2*nfreq*ndof
 
 
diff --git a/wecopttool/waves.py b/wecopttool/waves.py
index 7603f9b7..184048ab 100644
--- a/wecopttool/waves.py
+++ b/wecopttool/waves.py
@@ -52,8 +52,9 @@
 
 def elevation_fd(
     f1: float,
-    nfreq: int,
+    ifreq_end: int,
     directions: Union[float, ArrayLike],
+    ifreq_start: int = 1,
     amplitudes: Optional[ArrayLike] = None,
     phases: Optional[ArrayLike] = None,
     attr: Optional[Mapping] = None,
@@ -88,7 +89,7 @@ def elevation_fd(
     """
     directions = np.atleast_1d(degrees_to_radians(directions, sort=False))
     ndirections = len(directions)
-    freq = frequency(f1, nfreq, False)
+    freq = frequency(f1, ifreq_end)[ifreq_start:]
     omega = freq*2*np.pi
 
     dims = ('omega', 'wave_direction')
@@ -99,6 +100,7 @@ def elevation_fd(
               'freq': (dims[0], freq, freq_attr),
               'wave_direction': (dims[1], directions, dir_attr)}
 
+    nfreq = (ifreq_end - ifreq_start + 1)
     if amplitudes is None:
         amplitudes = np.zeros([nfreq, ndirections])
 
@@ -120,11 +122,12 @@ def elevation_fd(
 
 def regular_wave(
     f1: float,
-    nfreq: int,
+    ifreq_end: int,
     freq: float,
     amplitude: float,
     phase: Optional[float] = None,
     direction: float = 0.0,
+    ifreq_start: int = 1,
 ) -> DataArray:
     """Create the dataset for a regular wave.
 
@@ -149,7 +152,8 @@ def regular_wave(
 
     # attributes & index
     omega = freq*2*np.pi
-    tmp_waves = elevation_fd(f1, nfreq, direction)
+    nfreq = (ifreq_end - ifreq_start + 1)
+    tmp_waves = elevation_fd(f1, ifreq_end, direction, ifreq_start)
     iomega = tmp_waves.sel(omega=omega, method='nearest').omega.values
     ifreq = iomega/(2*np.pi)
 
@@ -176,7 +180,7 @@ def regular_wave(
 
     # wave elevation
     tmp = np.zeros([nfreq, 1])
-    waves = elevation_fd(f1, nfreq, direction, tmp, tmp, attrs)
+    waves = elevation_fd(f1, ifreq_end, direction, ifreq_start, tmp, tmp, attrs)
     waves.loc[{'omega': iomega}] = amplitude  * np.exp(1j*rphase)
 
     return waves