2213 lines
1.2 MiB
Plaintext
2213 lines
1.2 MiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "91e4c90f",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"Step 1: Clean out non-model simulation points: gap heights below 5mm"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"id": "e4ff9106",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import pandas as pd\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"import os\n",
|
||
|
|
"plt.rcParams['text.usetex'] = True"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 4,
|
||
|
|
"id": "32946653",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"<class 'pandas.core.frame.DataFrame'>\n",
|
||
|
|
"RangeIndex: 4471 entries, 0 to 4470\n",
|
||
|
|
"Data columns (total 6 columns):\n",
|
||
|
|
" # Column Non-Null Count Dtype \n",
|
||
|
|
"--- ------ -------------- ----- \n",
|
||
|
|
" 0 currL [A] 4471 non-null int64 \n",
|
||
|
|
" 1 currR [A] 4471 non-null int64 \n",
|
||
|
|
" 2 rollDeg [deg] 4471 non-null float64\n",
|
||
|
|
" 3 GapHeight [mm] 4471 non-null float64\n",
|
||
|
|
" 4 YokeForce.Force_z [newton] 4471 non-null float64\n",
|
||
|
|
" 5 YokeTorque.Torque [mNewtonMeter] 4471 non-null float64\n",
|
||
|
|
"dtypes: float64(4), int64(2)\n",
|
||
|
|
"memory usage: 209.7 KB\n",
|
||
|
|
"\n",
|
||
|
|
"After adding mirrored data:\n",
|
||
|
|
"<class 'pandas.core.frame.DataFrame'>\n",
|
||
|
|
"RangeIndex: 8281 entries, 0 to 8280\n",
|
||
|
|
"Data columns (total 6 columns):\n",
|
||
|
|
" # Column Non-Null Count Dtype \n",
|
||
|
|
"--- ------ -------------- ----- \n",
|
||
|
|
" 0 currL [A] 8281 non-null int64 \n",
|
||
|
|
" 1 currR [A] 8281 non-null int64 \n",
|
||
|
|
" 2 rollDeg [deg] 8281 non-null float64\n",
|
||
|
|
" 3 GapHeight [mm] 8281 non-null float64\n",
|
||
|
|
" 4 YokeForce.Force_z [newton] 8281 non-null float64\n",
|
||
|
|
" 5 YokeTorque.Torque [mNewtonMeter] 8281 non-null float64\n",
|
||
|
|
"dtypes: float64(4), int64(2)\n",
|
||
|
|
"memory usage: 388.3 KB\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"magDf = pd.read_csv(\"Ansys Results 12-9.csv\")\n",
|
||
|
|
"magDf.info()\n",
|
||
|
|
"\n",
|
||
|
|
"condition = magDf['GapHeight [mm]'] < 5\n",
|
||
|
|
"\n",
|
||
|
|
"magDf = magDf[~condition]\n",
|
||
|
|
"\n",
|
||
|
|
"nzRoll = magDf[magDf['rollDeg [deg]'] != 0]\n",
|
||
|
|
"\n",
|
||
|
|
"mirrored_data = nzRoll.copy()\n",
|
||
|
|
"mirrored_data['rollDeg [deg]'] = -nzRoll['rollDeg [deg]']\n",
|
||
|
|
"mirrored_data['currL [A]'], mirrored_data['currR [A]'] = nzRoll['currR [A]'].values, nzRoll['currL [A]'].values\n",
|
||
|
|
"mirrored_data['YokeTorque.Torque [mNewtonMeter]'] = -nzRoll['YokeTorque.Torque [mNewtonMeter]']\n",
|
||
|
|
"\n",
|
||
|
|
"magDf = pd.concat([magDf, mirrored_data], ignore_index=True)\n",
|
||
|
|
"\n",
|
||
|
|
"magDf = magDf.sort_values(['rollDeg [deg]', 'currL [A]', 'currR [A]', 'GapHeight [mm]']).reset_index(drop=True)\n",
|
||
|
|
"\n",
|
||
|
|
"print()\n",
|
||
|
|
"print(\"After adding mirrored data:\")\n",
|
||
|
|
"magDf.info()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "ecd1f124",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### After removing non-model rows, we have 4459 rows * ~2 for the mirrored roll angles, which matches the number of simulations conducted in Ansys."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "1a76017e",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAHCCAYAAAAXajikAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATS9JREFUeJzt3Xl8E3XeB/BPWqAg0E4DchfolCIgCCQprK7HalMPXFeFtOh6rBet7D67rldDPdbH9Sitx7rus0qKt+5BExDXExtQ12ulTSiC3JlyyE2TtOXomXn+CMm20JY2TTo5Pu/Xi9eSyWT6zWxJPv7mN7+vSpZlGUREREQRLk7pAoiIiIiCgaGGiIiIogJDDREREUUFhhoiIiKKCgw1REREFBUYaoiIiCgqMNQQERFRVGCoISIioqjAUENEQVFcXAyLxaJ0GQGLxPrz8vKgUqkgSVKn24hihYorChPFFrfbjdLSUpjNZpSVlbV5zmq1wmQyISsrC6IooqysDBkZGTAYDGc8ZnJyMkwmE3Jzc0NZftjVX1xcDEEQ/MfJz8/3P1dSUgKHwwEAyMrKgl6v79axz8RXd+uP8fa2EcWKPkoXQES9x263o6KiAm63G06n87Tn3W43rFYrLBYLRFGE0Wg8YyAAvF/eGo3G/wUeKuFWf3FxMQD4g5DVakVeXh5MJpN/pKSoqAgAYDQagx5qrFbracdsbxtRrGCoIYohGo0GGo2m08ssVVVV/pGHrvB9iZaXl4f8kke41V9YWIiqqir/Y71ej6ysLH+o0el0/ufS0tIgSRJEUezWz+hMWVkZsrKyzriNKFZwTg0R9YjdbodGo4FarYbb7Va6nG4LtH5JkuB2u9sNUL6gZDKZ4Ha7IUkSbDZbm0BjsVhgNBohSRJKSkpgsViQl5d32rGKi4tRUlLi/9PezznTNqJYwZEaImqjtLQUarUaTqcTDofDf/mkPSUlJf5LL2lpabBarb1VZod6q/6ORnUEQfCHo6KiIpSWlgIATCaTfx9fkJIkCdnZ2Vi9ejUkSTrt5/tGfXxhSKvVQqfTQaPR+MOSRqPx79/eNqJYwlBDRH6+L0Pfl2hJSQmys7NhNptP29ftdkOtVvtHKgRBaHeeS28Kh/p9gcp3zPYmHjudTv8cnvnz50MQhNPm9BiNRmg0mjajOzqdDlarFRqNhqM0RO1gqCGKUBaLBcuWLTvjfgUFBV3+L/dT53vk5OQgLy+v3cssRqMRaWlp/smy5eXl3bp8E+n1d6QrwcgXPHwTi9tTXFwMm83WZpskSUhLSwPA+TRE7ZKJKOaYzWZZo9G0u/1UAGSbzdZmm81ma3cbANnlcgW11vaEQ/0Oh0Nu7yMUgFxWVnbG17tcrnZf76tFEIR2j+2rWxTF095De9uIYgknChMRAO/lmOzs7DZzRXwjF6eOgPgugbTm20epRd96u35RFCEIQrv7d+USUEVFRYf7OZ1OqNXqNtssFov/7q/Wc2fsdjsAtLuNKNYw1BDFoPYukQiCgPz8/DYBoKSkBAaDoc2lm+Li4nbXfvHt0/rYvjt7gi1c6i8oKGgzuddisXR58b6ysrIOL6vp9frT3qPJZPLPDfLNyQG84aijbUSxhnNqiGKIJEn+uSx2ux1Go7HNirsFBQX+OSYAUF1d7f8ilSQJRqPRv0ZM65VzJUny32Xk+1+9Xg+r1Qqj0YicnJxurR0TKfXn5+e3aa9QXl7e5i6nM72XgoKCDp83m80oLi6GKIqQJAlms9lfgyiK0Ol0sFgs/tGe9rYRxRq2SSCikPJdngnmonO9KdLrJ4olvPxERCEV7FV0e1uk108USxhqiCikInGV4dYivX6iWMJQQ0Qh43a7I3qUI9LrJ4o1nFNDREREUYEjNURERBQVGGqIiIgoKsTUOjUejwf79u3D4MGDoVKplC6HiIiIukCWZdTV1WHUqFGIi+t4PCamQs2+ffuQkpKidBlEREQUgD179mDMmDEdPh9ToWbw4MEAvCclMTFR4WqIiIioK2pra5GSkuL/Hu9ITIUa3yWnxMREhhoiIqIIc6apI5woTERERFGBoYaIiIiiAkMNERERRQWGGiIiIooKDDVEREQUFRhqiIiIKCow1BAREVFUYKghIiKiqMBQQ0RERFGBoYaIiIiiAkMNERERRQWGGiIiIooKDDVERETUYzUnmvDWtzshy7JiNcRUl24iIiIKvkN19fjVa+XYvL8WJxpbkHdJmiJ1MNQQERFRwPY4j+PmV7/DrurjGDooARdPPFuxWhhqiIiIKCDbDtbhlle/w8HaBoxJHoB37pyN8UMHKlYPQw0RERF127rdLtz+Rjncx5swcfggvHXHbIxI6q9oTQw1RERE1C1fbT+C3LcrcLyxBTNSBLxxewaEs/opXRZDDREREXXdJxv343f/qERjiwcXThgK0y1aDEwIjzgRHlUQERFR2Cst34NFK76HRwaumjoCL9wwAwl94pUuy4+hhoiIiM6o5N8OPP3RFgDAfF0Knp47DfFxKoWraouhhoiIiDokyzKeWbUVL33uAADkXSJi0ZWToFKFV6ABGGqIiIioAy0eGY++txF//243AMB45SQs/JkyC+t1BUMNERERnaax2YN7Syvx4ff7oVIBT103Db+cPVbpsjrFUENERERtHG9sxt3v2PHvbYfRN16FP82fgZ+fN0rpss5IsVBjt9thtVoBAOXl5Vi6dCkEQfA/BwAajQaSJMHtdkOj0QAAJEmCxWKBKIqQJAm5ubn+1xEREVHP1Bxvwh1vlsO2y4UBfeOx5BYtLlGw9UF3KBZqrFYr8vPzAQDFxcXIzMyEzWYDAJhMJpSUlAAA9Ho9zGaz/3XZ2dn+/SRJwoIFC9o8T0RERIE5VFuPW19biy0H6pDYvw9ev30WtOOSlS6ry+KU+KF2ux2FhYX+xwaDAXa7HZIkAQC0Wi1cLhdcLhfKysr8IzG+531EUfSP9hAREVHgdlcfh2HJt9hyoA5nD05A6d3nR1SgARQaqdFoNFi6dKn/sdvtBgCo1Wr/tvYuKVmt1jb7+F5jt9v9l6daa2hoQENDg/9xbW1tDysnIiKKPlsPeBtTHqprwFj1WXjnztkYO+QspcvqNsUuPxkMBv/fly1bBr1e7w8ybrcbFosFgHe+TV5eHkRR9IefUzmdzna3FxYW4vHHHw9q3URERNHEvtuF218vR82JJpwzfDDevnMWhiUq25gyUIrf/eQLML55MgDaTP4VRRFZWVlwOBydHqM9BQUFuO+++/yPa2trkZKSEpS6iYiIIt2X2w8j9y0bTjS1YOZYAa/fFh6NKQOleKgxGo1t5s0A3rkzvstJvrucJEmCIAinjco4nc4O735KSEhAQkJCqEonIiKKWB9t2I97/rkOTS0yLkr3NqY8q5/isaBHFJko7FNcXAyj0ei/tOR2u2G325GZmXnavmq1Gnq9vt3j6HS6UJdKREQUNf65djf+5+92NLXIuHraSLzyK13EBxpAwVBjsVig0Wj8gaa0tBSCIEAURRQVFfn3s1qtMBgM/udakyQJOp2O69QQERF10ZIvHFi0YgM8MnDjrBS8eOPMsOq03RMqWZbl3v6hkiQhLa1t7whBEOByuQD8d2E+QRDgcDjahBxJkmAymZCRkYHy8nIUFBR0OdTU1tYiKSkJNTU1SExMDNr7ISIiCneyLKPok61Y8oV3jurCn6Uh/4pzwrIx5am6+v2tSKhRCkMNERHFohaPjEdWbsA/1u4BABRcNQl5l4RvY8pTdfX7O/IvoBEREVGHGppbcN+y9fhww37EqYCnr5+GG2aFd2PKQDHUEBERRaljDc24+x0bvtx+BH3jVfjzDTMxZ9pIpcsKGYYaIiKiKOQ+3ojb3yjHut1unNUvHqZbtLgoPTIaUwaKoYaIiCjKHKqtxy2vrsXWg3VIGtAXr9+eAc3YyOrjFAiGGiIioiiyq/oYbn71O+xxnsCwwQl4+87ZOGfEYKXL6hUMNURERFFiy4Fa3PLqWhyua8C4Id7GlCnqyGtMGSiGGiIioihg2+XE7a+Xo7a+GZNGDMZbd87CsMGR2ZgyUAw1RER
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 640x480 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"first_subset = magDf.iloc[0:13]\n",
|
||
|
|
"x = first_subset[\"GapHeight [mm]\"]\n",
|
||
|
|
"x = 1/x\n",
|
||
|
|
"y = first_subset[\"YokeForce.Force_z [newton]\"]\n",
|
||
|
|
"plt.plot(x, y)\n",
|
||
|
|
"plt.title(r\"$-15A, -15A, 0^\\circ roll$\")\n",
|
||
|
|
"plt.xlabel(\"Inverse Gap Height (1/mm)\")\n",
|
||
|
|
"plt.ylabel(\"Force on magnetic yoke (N)\")\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"# Experimented here with ways to represent gap-height-to-force relationship."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 34,
|
||
|
|
"id": "1ba9a1b6",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"<class 'pandas.core.frame.DataFrame'>\n",
|
||
|
|
"RangeIndex: 8281 entries, 0 to 8280\n",
|
||
|
|
"Data columns (total 7 columns):\n",
|
||
|
|
" # Column Non-Null Count Dtype \n",
|
||
|
|
"--- ------ -------------- ----- \n",
|
||
|
|
" 0 currL [A] 8281 non-null int64 \n",
|
||
|
|
" 1 currR [A] 8281 non-null int64 \n",
|
||
|
|
" 2 rollDeg [deg] 8281 non-null float64\n",
|
||
|
|
" 3 GapHeight [mm] 8281 non-null float64\n",
|
||
|
|
" 4 YokeForce.Force_z [newton] 8281 non-null float64\n",
|
||
|
|
" 5 YokeTorque.Torque [mNewtonMeter] 8281 non-null float64\n",
|
||
|
|
" 6 invGap 8281 non-null float64\n",
|
||
|
|
"dtypes: float64(5), int64(2)\n",
|
||
|
|
"memory usage: 453.0 KB\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/html": [
|
||
|
|
"<div>\n",
|
||
|
|
"<style scoped>\n",
|
||
|
|
" .dataframe tbody tr th:only-of-type {\n",
|
||
|
|
" vertical-align: middle;\n",
|
||
|
|
" }\n",
|
||
|
|
"\n",
|
||
|
|
" .dataframe tbody tr th {\n",
|
||
|
|
" vertical-align: top;\n",
|
||
|
|
" }\n",
|
||
|
|
"\n",
|
||
|
|
" .dataframe thead th {\n",
|
||
|
|
" text-align: right;\n",
|
||
|
|
" }\n",
|
||
|
|
"</style>\n",
|
||
|
|
"<table border=\"1\" class=\"dataframe\">\n",
|
||
|
|
" <thead>\n",
|
||
|
|
" <tr style=\"text-align: right;\">\n",
|
||
|
|
" <th></th>\n",
|
||
|
|
" <th>currL [A]</th>\n",
|
||
|
|
" <th>currR [A]</th>\n",
|
||
|
|
" <th>rollDeg [deg]</th>\n",
|
||
|
|
" <th>GapHeight [mm]</th>\n",
|
||
|
|
" <th>YokeForce.Force_z [newton]</th>\n",
|
||
|
|
" <th>YokeTorque.Torque [mNewtonMeter]</th>\n",
|
||
|
|
" <th>invGap</th>\n",
|
||
|
|
" </tr>\n",
|
||
|
|
" </thead>\n",
|
||
|
|
" <tbody>\n",
|
||
|
|
" <tr>\n",
|
||
|
|
" <th>0</th>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-4.0</td>\n",
|
||
|
|
" <td>6.0</td>\n",
|
||
|
|
" <td>260.518631</td>\n",
|
||
|
|
" <td>-6976.851677</td>\n",
|
||
|
|
" <td>0.166667</td>\n",
|
||
|
|
" </tr>\n",
|
||
|
|
" <tr>\n",
|
||
|
|
" <th>1</th>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-4.0</td>\n",
|
||
|
|
" <td>8.0</td>\n",
|
||
|
|
" <td>163.529872</td>\n",
|
||
|
|
" <td>-3335.704070</td>\n",
|
||
|
|
" <td>0.125000</td>\n",
|
||
|
|
" </tr>\n",
|
||
|
|
" <tr>\n",
|
||
|
|
" <th>2</th>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-4.0</td>\n",
|
||
|
|
" <td>9.0</td>\n",
|
||
|
|
" <td>136.554807</td>\n",
|
||
|
|
" <td>-2506.094902</td>\n",
|
||
|
|
" <td>0.111111</td>\n",
|
||
|
|
" </tr>\n",
|
||
|
|
" <tr>\n",
|
||
|
|
" <th>3</th>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-4.0</td>\n",
|
||
|
|
" <td>10.0</td>\n",
|
||
|
|
" <td>117.403213</td>\n",
|
||
|
|
" <td>-1959.725693</td>\n",
|
||
|
|
" <td>0.100000</td>\n",
|
||
|
|
" </tr>\n",
|
||
|
|
" <tr>\n",
|
||
|
|
" <th>4</th>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-15</td>\n",
|
||
|
|
" <td>-4.0</td>\n",
|
||
|
|
" <td>10.5</td>\n",
|
||
|
|
" <td>109.107025</td>\n",
|
||
|
|
" <td>-1759.467304</td>\n",
|
||
|
|
" <td>0.095238</td>\n",
|
||
|
|
" </tr>\n",
|
||
|
|
" </tbody>\n",
|
||
|
|
"</table>\n",
|
||
|
|
"</div>"
|
||
|
|
],
|
||
|
|
"text/plain": [
|
||
|
|
" currL [A] currR [A] rollDeg [deg] GapHeight [mm] \\\n",
|
||
|
|
"0 -15 -15 -4.0 6.0 \n",
|
||
|
|
"1 -15 -15 -4.0 8.0 \n",
|
||
|
|
"2 -15 -15 -4.0 9.0 \n",
|
||
|
|
"3 -15 -15 -4.0 10.0 \n",
|
||
|
|
"4 -15 -15 -4.0 10.5 \n",
|
||
|
|
"\n",
|
||
|
|
" YokeForce.Force_z [newton] YokeTorque.Torque [mNewtonMeter] invGap \n",
|
||
|
|
"0 260.518631 -6976.851677 0.166667 \n",
|
||
|
|
"1 163.529872 -3335.704070 0.125000 \n",
|
||
|
|
"2 136.554807 -2506.094902 0.111111 \n",
|
||
|
|
"3 117.403213 -1959.725693 0.100000 \n",
|
||
|
|
"4 109.107025 -1759.467304 0.095238 "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 34,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# sooooo it looks like we can invert the entire gap height column and see how that works out...\n",
|
||
|
|
"magDf[\"invGap\"] = magDf[\"GapHeight [mm]\"].transform(lambda x: 1/x)\n",
|
||
|
|
"magDf.info()\n",
|
||
|
|
"magDf.head()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "9cd13fd5",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Let's try fitting a polynomial to this..."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "bd16a120",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Force Coeffs: [-8.65013483e-02 2.83478682e-01 2.86948944e-01 -2.40620478e-08\n",
|
||
|
|
" 3.53295625e+02 -1.79173973e-03 -1.72965671e-03 -2.24271795e-02\n",
|
||
|
|
" -1.14467614e+01 -2.16073374e-03 2.24271792e-02 -1.16266321e+01\n",
|
||
|
|
" 3.17540506e-01 4.74016321e-07 6.56934801e+03 -4.22077144e-05\n",
|
||
|
|
" -1.06779010e-05 -2.15844424e-04 1.18908944e-01 -7.64610873e-06\n",
|
||
|
|
" 5.14110976e-12 1.01686168e-01 -1.01983866e-02 1.54885341e-01\n",
|
||
|
|
" -1.03778353e+02 -4.02599676e-05 2.15844412e-04 1.29958897e-01\n",
|
||
|
|
" -1.01471001e-02 -1.54885342e-01 -1.00615042e+02 -4.31582506e-10\n",
|
||
|
|
" 2.06864028e+00 -5.14832777e-06 -4.19530520e+04 -3.84317353e-07\n",
|
||
|
|
" 3.91430632e-07 -7.41836487e-07 7.93047251e-04 -4.85823477e-07\n",
|
||
|
|
" -2.53747638e-06 1.17129705e-04 -2.29814329e-04 4.03023222e-04\n",
|
||
|
|
" 2.97296739e-01 -1.05011736e-07 2.53747734e-06 -3.78155618e-06\n",
|
||
|
|
" -4.82390183e-05 2.17240115e-12 -5.68285028e-01 3.49652512e-03\n",
|
||
|
|
" 4.38290920e-01 -8.84140414e+00 7.71042932e+02 1.32237943e-07\n",
|
||
|
|
" 7.41839504e-07 7.27510024e-04 -2.31026411e-04 -4.03023131e-04\n",
|
||
|
|
" 1.91281517e-01 -3.49652508e-03 4.35555849e-01 8.84140414e+00\n",
|
||
|
|
" 7.49444072e+02 -3.20760957e-02 -7.80704790e-10 -1.33481028e+02\n",
|
||
|
|
" 3.32901902e-05 7.71176022e+04 3.63954200e-08 1.21055428e-08\n",
|
||
|
|
" 1.38326925e-07 -2.26843332e-05 2.02887840e-08 -1.26534403e-07\n",
|
||
|
|
" -5.30119171e-06 8.63672085e-08 4.12729183e-05 -8.00925222e-04\n",
|
||
|
|
" 2.31832331e-08 2.02504680e-13 1.06948464e-05 -7.75304869e-07\n",
|
||
|
|
" 3.44199510e-05 7.75195644e-04 9.67991899e-06 3.94281516e-03\n",
|
||
|
|
" 1.52654027e-01 -5.85408049e-01 2.66528575e-08 1.26534403e-07\n",
|
||
|
|
" -1.16260075e-07 -8.31899683e-07 -3.44199134e-05 1.37334995e-03\n",
|
||
|
|
" -7.83373366e-13 7.56972747e-04 3.87534095e-10 1.67361829e+00\n",
|
||
|
|
" -1.53624144e-04 -5.22167175e-02 -4.26268983e+00 -1.86825486e+01\n",
|
||
|
|
" -1.80981191e+03 3.91171451e-08 -1.38326669e-07 -2.90007777e-05\n",
|
||
|
|
" 5.72601948e-08 -4.12729382e-05 -4.53297467e-04 -9.67991778e-06\n",
|
||
|
|
" 3.95070566e-03 -1.52654026e-01 -2.35783054e-01 -1.49754210e-04\n",
|
||
|
|
" 5.22167174e-02 -4.24980327e+00 1.86825486e+01 -1.75997408e+03\n",
|
||
|
|
" 3.41748851e-12 4.47231699e-01 1.88369113e-09 1.08183457e+03\n",
|
||
|
|
" -8.46779848e-05 4.14731651e+04]\n",
|
||
|
|
"Torque Coeffs: [ 8.37786404e-01 8.73852020e+00 -8.81670143e+00 -6.05820755e+01\n",
|
||
|
|
" -9.82722151e+01 -8.14346739e-02 -8.55581205e-03 -1.26586317e+00\n",
|
||
|
|
" -2.28782708e+02 7.92025577e-02 -1.26586317e+00 2.32329851e+02\n",
|
||
|
|
" 1.59724819e-01 2.40322138e+02 1.91664184e+03 -1.60379946e-03\n",
|
||
|
|
" 1.33688429e-03 -1.09883437e-02 4.99001714e+00 -1.13086428e-03\n",
|
||
|
|
" -5.01219559e-04 2.95246925e-01 -5.53027170e-01 1.40436065e+01\n",
|
||
|
|
" -7.51330187e+03 1.58512910e-03 -1.09883435e-02 -5.01157061e+00\n",
|
||
|
|
" 5.39147403e-01 1.40436064e+01 7.47720714e+03 8.04333337e+00\n",
|
||
|
|
" -9.27394310e-01 5.61703432e+04 -1.36591792e+04 -8.85041568e-05\n",
|
||
|
|
" 5.78299126e-05 -3.32035777e-04 1.42291379e-02 9.30848750e-06\n",
|
||
|
|
" -1.81942915e-04 4.52448003e-03 -1.03937242e-02 -1.18018293e-01\n",
|
||
|
|
" 1.97297590e+01 -5.50802583e-05 -1.81942945e-04 -3.44846992e-03\n",
|
||
|
|
" 5.10224715e-05 -1.04187796e-01 -3.02785121e+00 2.03017250e-01\n",
|
||
|
|
" 2.26755131e+01 -3.31637758e+02 4.95131780e+04 9.86642026e-05\n",
|
||
|
|
" -3.32035797e-04 -1.38739075e-02 1.04393208e-02 -1.18018293e-01\n",
|
||
|
|
" -1.90316038e+01 2.03017250e-01 -2.25845770e+01 -3.31637758e+02\n",
|
||
|
|
" -4.94386582e+04 -4.28915364e-03 -2.85648016e+02 -2.71180777e+00\n",
|
||
|
|
" -4.00600715e+05 2.87296873e+04 1.94952853e-06 -5.18824865e-06\n",
|
||
|
|
" -7.73369810e-06 5.36212900e-04 1.10846270e-06 -1.84292257e-06\n",
|
||
|
|
" -4.65905674e-04 6.42855031e-06 3.08079390e-03 -1.09734095e-03\n",
|
||
|
|
" -1.62663865e-06 6.42939658e-06 -1.00143627e-04 -2.50122237e-06\n",
|
||
|
|
" 1.71620024e-03 -2.14975506e-02 7.20307314e-04 1.76411710e-01\n",
|
||
|
|
" 8.07086586e+00 -3.68279950e+01 4.55521791e-06 -1.84292611e-06\n",
|
||
|
|
" 4.37141266e-04 -1.60125131e-06 1.71620050e-03 1.99693938e-02\n",
|
||
|
|
" 1.48517227e-04 -1.19886653e-03 1.78969383e+00 9.92388924e+00\n",
|
||
|
|
" -7.19412178e-03 -2.98153181e+00 -2.04344115e+02 -1.80247522e+03\n",
|
||
|
|
" -1.11849627e+05 -1.66014902e-06 -7.73369800e-06 -7.23358309e-04\n",
|
||
|
|
" -4.16903359e-06 3.08079402e-03 -6.84048059e-03 7.20307304e-04\n",
|
||
|
|
" -1.75664664e-01 8.07086586e+00 3.49676892e+01 7.76959664e-03\n",
|
||
|
|
" -2.98153181e+00 2.03737293e+02 -1.80247522e+03 1.12167507e+05\n",
|
||
|
|
" 4.69772092e-02 -2.00823262e-02 2.13341632e+03 5.37309061e+01\n",
|
||
|
|
" 1.24121750e+06 1.52681158e+04]\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"from sklearn.preprocessing import PolynomialFeatures\n",
|
||
|
|
"from sklearn.linear_model import LinearRegression\n",
|
||
|
|
"\n",
|
||
|
|
"# The following was written by AI - see [4]\n",
|
||
|
|
"# 1. Load your Ansys CSV - using invGap for modeling\n",
|
||
|
|
"X = magDf[['currL [A]', 'currR [A]', 'rollDeg [deg]', 'invGap']]\n",
|
||
|
|
"y = magDf[['YokeForce.Force_z [newton]', 'YokeTorque.Torque [mNewtonMeter]']]\n",
|
||
|
|
"\n",
|
||
|
|
"# 2. Create Features (e.g. z^2, z^3, z*IL, etc.)\n",
|
||
|
|
"poly = PolynomialFeatures(degree=5) \n",
|
||
|
|
"X_poly = poly.fit_transform(X)\n",
|
||
|
|
"\n",
|
||
|
|
"# 3. Fit\n",
|
||
|
|
"model = LinearRegression()\n",
|
||
|
|
"model.fit(X_poly, y)\n",
|
||
|
|
"\n",
|
||
|
|
"# 4. Extract Equation\n",
|
||
|
|
"print(\"Force Coeffs:\", model.coef_[0])\n",
|
||
|
|
"print(\"Torque Coeffs:\", model.coef_[1])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "dd6861d0",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Compare Model Predictions with Raw Data"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "dc79d9fe",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABW4AAAHpCAYAAAASzqVtAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAA5CpJREFUeJzs3XlcVXX+x/HXvaxXWS4g7pRetLJNBWwzKwPbm2kBrfm1Tglly8zUKDpNM9NMjUJT02qh1Uw1SwpNi+2i7atA0l7K1cQd4V4QvKz3/P4gjiKoIMjlwvv5ePB4eL/L5XPPwfr64Xu+H4thGAYiIiIiIiIiIiIi0mtYfR2AiIiIiIiIiIiIiLSmxK2IiIiIiIiIiIhIL6PErYiIiIiIiIiIiEgvo8StiIiIiIiIiIiISC+jxK2IiIiIiIiIiIhIL6PErYiIiIiIiIiIiEgvo8StiIiIiIiIiIiISC+jxK2IiIiIiIiIiIhIL6PErYjsV35+PhaLZZ9f8fHxvg6xS4qKitr9XImJiWRnZx/S791ybZ1Op9m2aNEi0tLSDun3PZD24uoNEhMTsVgsFBUV+TqUA8rOziYqKsrXYYiIiEg78vPziYqK6vBXb1sTHQpaE7elNXHXaU0s0nWBvg5ARPzDnDlzmDFjBgAVFRVme3R0tK9C6lY5OTlMnz4daP58+fn5ZGZmsmTJEgoLC3ssjsLCQvLy8nrs+/kLp9NJUVERdrudnJwccnJyfB2SiIiI+KmkpCRyc3NbteXm5rJo0SKysrJISEho1edwOHoyPJ/Smrh305pYpP/RjlsR6ZD4+HgSEhJISEggJSXF/Np7YeuvoqOjsdvt2O12HA4H6enpFBYWUlRURGZmZo/FkZOTg2EYnZ6Xl5eH2+3u/oB6iby8PFJSUpg+fTpLly7t8nv15WslIiIi+2e321utZ1NSUkhMTARo056SkuLjaHuW1sS9m9bEIv2PErciIvvgcDhITU0lPz/f16EcUFpaGgUFBb4O45DJyckhLS2NtLQ03G53l+5JX79WIiIiIt1Ja+LeQ2tikf5HiVsRkQPQb6J9q6ioCKfT2WrXy96PN4qIiIjIoaU1sW9pTSzSPylxKyLdpqioiGnTphEVFUV8fHybQgbZ2dlMmzYNgIyMDKKioswF4J5zo6KimDZtWqsD9/Py8syD+BMTE1v9dtnpdGKxWFi0aFG3fp6W32Lv+Yjc/j7D/mKE5iIL8fHxREVFkZaW1m6hg/YO8N/ftUlLS8NisQAwbdo0s5DEnrojrn3JzMxst+DA3sUc8vPzzRjau7/7s2TJEhwOh3m+XEpKyn7v9b6u14GuVXufpb2iFE6nk7S0NKKiosxr6g/FIUREROTgdWWdm5eX12qtlZeX12Z90dF1SMv79eS6WGviA9OaWGtikUNFiVsR6ZCMjIx2K822LNBaFiEJCQmsWLGCrKws81GePTmdTnOhNG/ePOx2O0VFRSQmJmK328nNzSU3N5eEhASWLFkC7K4qm5GRQWFhITNmzGDatGnmwiE6Opo5c+aQlJTULZ+1ZXHactZZVlbWAT/DgWJctGgRGRkZ5vWZMWNGh84JO9C1Wbx4sVkoIjc3F5fLhcvlMucfqrhazJgxo93HtHJyckhISMDhcOB2u5k2bRozZsygpKTE/Ax7Frnbn0WLFpGammq+bvmZaq9gxf6u14GuVUfl5eURHR1Nbm4uJSUlOBwOkpOTtQtFRESkj+rKOjcvL4+0tDQcDge5ublMmzaNmTNnHnQsPbku1ppYa+L90ZpYpIcYIiL7sXz5cgMwsrKyjJKSkjZfLRwOhzFnzpxWc0tKSgzAyM3NNQzDMLKysgzASElJaTXO4XAYqamp7X5/l8tlAEZOTk6r9jlz5hjp6eld/nyFhYUG0ObLbrcbqamphsvlajW+vc/QkRhb3m9Pubm5BtDqOmZlZRl2u918vb9rs/f3X758ebvt3RHX/jgcjjb3ouVnxjB2/wwdjJb7U1hYaLa1fK72rsuBrte+rpVhNF+XPa/9nrHv71q0d533vo8iIiLSu+Xk5LRZc7ToyjrXbrcbCQkJrdrmzJnTZn3RkXXIoVwXa02sNXELrYlFehftuBWRDmmpLLv3F+w+bykjI6PVHIfD0eo34S1ycnLMPzudznbntmg5MH/vHb/Z2dndeph+VlaW+dvmlq/c3Fzsdnu74/f8DAeK0el04na79/kZ9+VA1+ZADlVce0tNTW1V1bblt/7p6ekA5o6PxMREFi1a1KnHznJycsyfPbfbbf4GPyEhoc3ugq5er4PV8jNSUlLSo99XREREDr2urnPdbnebnaotRwx0Vk+si7UmPnhaE9sBrYlFupsStyLSZS2Loejo6DZ9DoejzaKkJeELmH17trXH5XJhGEarr5ZHfLqD3W5v9XUg7cW7rxg7+hn3drDzDnVce8vIyGj1aNiSJUtISUkxr6PdbjfvVUZGBvHx8UybNq1Dj1EtXboUt9ttns3V8tVyftaeC9Xu+jwdkZ+fT1pamnkOmoiIiPRNPbHO7axDuS7WmvjgaU2sNbHIoaDErYh0WcuioL3fGjudzlaLhr0XgPub25F+X+jsZ2hZ6Hf2M3T1sx+quNr7PgkJCWZV25az3PaUkJBAYWEhLpeLnJwcCgoKmD9//n7fNz8/H7fbTUlJSZtFdss5XHvu8uipn5W0tDTS0tKYNm0ay5cvP6gzwURERMQ/dGWd211rrY7E4gtaE7f9PloTi0h3U+JWRLqs5TfJey4YoPnRsqKiImbMmLHPuS1HLuw9F5oLIrQsgNpb0PSWg+8PFGNCQoJZGGBPBypEcKBr06Jl0bz39ThUcbVnxowZLF261Pxt//Tp09sdZ7fbSU9PJyUl5YBVZ1uKKLS3W8But5OSkmIuZKFj12tf12rvcS32jtHtdpOXl8fixYtJT0/vkZ0MIiIi4jtdWee2rLXaK+rVngOtQ3r7ulhrYq2JRaT7KXErIt1i8eLFZjXW/Px8Fi1aRHJyMikpKa2qn7YnJyeHvLw8c27Ln5OTk833bmkrKioiPz+fjIwMsyKv2+0mMzPzgIueQ+lAMc6bN49FixaZcbb8+UAOdG1a2O12lixZQlFRkfmY1qGMa2/p6em43W7mz59Pampqqx0YeXl5JCYmkpeXR1FREXl5eeTn5x/wfLelS5fu9x9DLed27XmWWEeu176uVXx8PLB7V0NeXl6bBX7LY4Pz588nPz+foqKiNjspREREpG/pyjp33rx55vqrZa3V3tmjHVmHtMTSm9fFWhNrTSwi3ezQ1z8TEX/WUkF07yqs7SksLDRSUlIMwHA4HGYF1Rb7qyraMtdut5tVXfesWlpSUtKqPz093axu21LVtyMxtvd9OzN3f59hfzHuObelLzc310hISGh3zN4x7u/atMxrue57f9/uiKsjWu59S3XlPb9/enq64XA4zBj3rsy8t5afuwPFQDvVmw90vfZ1rVwul5GQkGBWUE5PTzeWL19uJCQktJqfm5trvndCQoKRk5NjpKSkqIKuiIiIH8vJyTEAo7CwsN3+rqxzW/pa1h4t32vP9UVH1yGGcWjWxVoTa02sNbFI72QxDMPouTSxiIiIiIiISP/VcvZpSUmJHjEXEZH90lEJIiIiIiIiIiIiIr2MErciIiIiIiIiIiIivYwStyIiIiIiIiIiIiK9jM64FREREREREREREelltONWREREREREREREpJdR4lZERERERERERESklwn0dQC+4vV62bx5M+Hh4VgsFl+HIyIiItLvGYbBzp07GT58OFar9hd0B615RURERHqXzqx5+23idvPmzcTFxfk6DBERERHZS2lpKSNHjvR1GH2C1rwiIiIivVNH1rz9NnEbHh4ONF+kiIiIDs/zer2UlZURGxurnSB+RPfNP+m++SfdN/+k++af+tp9q6qqIi4uzlynSdcd7JrX3/S1vwu+pGvZvXQ9u4+uZffRtexeup7dp79cy86seftt4rblUbGIiIhOJ25ra2uJiIjo0z9EfY3um3/SffNPum/+SffNP/XV+6ZH+rvPwa55/U1f/bvgC7q
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1400x500 with 2 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Force R² Score: 0.999873\n",
|
||
|
|
"Torque R² Score: 0.999661\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"y_pred = model.predict(X_poly)\n",
|
||
|
|
"\n",
|
||
|
|
"force_pred = y_pred[:, 0]\n",
|
||
|
|
"torque_pred = y_pred[:, 1]\n",
|
||
|
|
"\n",
|
||
|
|
"force_actual = y.iloc[:, 0].values\n",
|
||
|
|
"torque_actual = y.iloc[:, 1].values\n",
|
||
|
|
"\n",
|
||
|
|
"fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
|
||
|
|
"\n",
|
||
|
|
"# Force comparison\n",
|
||
|
|
"axes[0].scatter(force_actual, force_pred, alpha=0.5, s=10)\n",
|
||
|
|
"axes[0].plot([force_actual.min(), force_actual.max()], \n",
|
||
|
|
" [force_actual.min(), force_actual.max()], \n",
|
||
|
|
" 'r--', linewidth=2, label='Perfect Fit')\n",
|
||
|
|
"axes[0].set_xlabel('Actual Force (N)', fontsize=12)\n",
|
||
|
|
"axes[0].set_ylabel('Predicted Force (N)', fontsize=12)\n",
|
||
|
|
"axes[0].set_title('Force: Predicted vs Actual', fontsize=14)\n",
|
||
|
|
"axes[0].legend()\n",
|
||
|
|
"axes[0].grid(True, alpha=0.3)\n",
|
||
|
|
"\n",
|
||
|
|
"# Torque comparison\n",
|
||
|
|
"axes[1].scatter(torque_actual, torque_pred, alpha=0.5, s=10)\n",
|
||
|
|
"axes[1].plot([torque_actual.min(), torque_actual.max()], \n",
|
||
|
|
" [torque_actual.min(), torque_actual.max()], \n",
|
||
|
|
" 'r--', linewidth=2, label='Perfect Fit')\n",
|
||
|
|
"axes[1].set_xlabel('Actual Torque (mN·m)', fontsize=12)\n",
|
||
|
|
"axes[1].set_ylabel('Predicted Torque (mN·m)', fontsize=12)\n",
|
||
|
|
"axes[1].set_title('Torque: Predicted vs Actual', fontsize=14)\n",
|
||
|
|
"axes[1].legend()\n",
|
||
|
|
"axes[1].grid(True, alpha=0.3)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"from sklearn.metrics import r2_score\n",
|
||
|
|
"print(f\"Force R² Score: {r2_score(force_actual, force_pred):.6f}\")\n",
|
||
|
|
"print(f\"Torque R² Score: {r2_score(torque_actual, torque_pred):.6f}\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "5a5d0634",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Visualize a specific subset (like the first 13 rows)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "837a814f",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA08AAAIjCAYAAADbfyCPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjrZJREFUeJzs3XlcXPW9//HXmWFPgGEL2UjCkLgl0TCAxmjqkklcu9hAYhfb3mpArV2urcH09i7e2zZCbXvv79pWEq1tbW8NoG1dayBqjUuUJVGTqDEM2XdgWJKwzvn9gYwhQAIZYAZ4Px8PhDnnO+d8znwnOB++3/P5GqZpmoiIiIiIiMgZWfwdgIiIiIiIyEig5ElERERERKQflDyJiIiIiIj0g5InERERERGRflDyJCIiIiIi0g9KnkRERERERPpByZOIiIiIiEg/KHkSERERERHpByVPIiKDzO12s3jxYn+HMaRG6zXm5+cTExPT52NfuFwuDMPAMAwqKyvP2DYlJQXDMMjNzR2UcwOUlpZiGAYul6vfzxnM6x8saWlp/XoNRUSGgpInERmT1qxZ0+cH07S0NIqLi3G73bjdbvLz88nPz+/3sVesWEF5eflghXrOAvkaXS4XOTk53vOuWbOm1/1dXwP5wD8SrFu3rs99lZWVo+56B4vL5aKyshKbzUZBQYG/wxGRMSjI3wGIiAwXl8tFXl4eAIWFhWRnZ/farrKykqysLO/j7Ozsfn9QKy0tpbKyErfb7XO852IkXKPL5SItLY3q6mpsNhsAubm55Ofns3LlSu/joqIi73OysrK6PR7JnE4na9as8fbT6datW4fT6aS0tHSYIwt8xcXFOJ1O7HY7hYWFSqBEZNhp5ElExgy73U5BQQEFBQXY7fY+23UlEgUFBVRVVQ3oA1pJSYl3tMcfCdRIuMa8vDyys7O9iRPAqlWrvMd0uVw9pgRmZGT4LSEdbFlZWbjd7j6TozVr1nRLbOVTBQUFZGVlnfU1FBEZKkqeREROk5KSQnZ2NtnZ2WdMQE63Zs0acnJyiI2NBQjoqVf+vMbCwkJSUlK6betKpEpLS4mNjaWqqqrb/pqamm7J1khmt9txOBy9jqSVlpbidrtZtmyZHyILbF3TGZ1OJ06nE2DUjEaKyMih5ElEpBddf9Xu703pbrebqqoq7Ha7Nxmpra0dyhB95o9r7LrHqreEzWazee9nAbyjCsXFxQM6R35+vnfkKicnh5iYGO+oVWVlJYsXLyYmJoaUlJQB3efVpavww+n3aQ1ETk5Or88vKirC6XT2mSgOJP41a9aQkpJCTEwMWVlZfSa6xcXF3iIMaWlp5zSak5ub22thidOLVJSWlnrPFRMTw+LFi/v9/lu3bl23917X9Mfe5Ofnk5aW5n29uq7t9HP1Fc9wXI+IjExKnkRETlNSUkJpaSnp6ekA/fpAtHr1alatWgV8OooSyCNP/rrGM7WPjY2lpqYG6Jza53a7vVP5+ro/6Ezn6UoEVq1ahc1m837QdTgcbNiwgby8PO80sIGIjY1l5cqV3tfuXHSNLJ2eGJ5pyt5A4u8aIexqu3z58l6Lh3SdLycnh4qKCpYvX87ixYsH3K/Lly/vdRpdQUEBDocDu93urdC4fPlyqqqqKCoqwuFw9DsBX7NmDZmZmd7HXdfdW3JdU1NDZWUlK1asIDc3l5KSEtxuN4sWLfK2OVM8OTk5uN3uHscezOsRkRHKFBEZgxwOh7ly5cpe91VVVXV7XFRUZNrt9j6PVVVVZRYUFHTbBph5eXm+B+qDQLzGiooKEzBLSkp67LPb7WZ2dvaAjtebvLw8EzCdTmeP45/+elRVVZmAWVRU5H2uzWbrdqxTH/ui61xd1+50Os3MzEzv/qKiIhMw6+rqTNPsfH1Pjbc/8Xex2Wzdjn3q8bv6vq6uzgR69OvKlSu9/TCQ6++t/059j5SUlJjn+rGj631TUVHh3dYV/+nX2XUNp7/PCgoKur2+Z4vH6XSaDoejx/V0vV6+XI+IjFwaeRIROc3pU8ocDgcul6vP6UxdBRBOd/p9O4HEX9d4pvuWBvsv9qcWwei6XyYnJ6dbm677j85UOnyo5OTkdBvZWLduHZmZmb2+RgOJ3+Vy4Xa7e7Q9XVep+ZycHO/6U4ZhkJ+ff05l6DMzMyksLPQ+7rq2rvdN10hdWloaa9asGdDoVkFBATabzTvi0zUN0+FwnHFaZ9e9UUCPqaZniycnJ6db2fiuKYJdo4a+XI+IjFwqVS4iI0rXVJz+Wrt2LQ6Ho9/tc3NzWb58ebfnnKk4QtcHt9OnRNlstnP+MDWar7HrPL1VznO73YNaFOLUBLErGeg6/+nt/PHBt2sKWnFxMZmZmRQXF/dZAGEg8Xf93N9CIHV1dYPyunet21VaWorT6fSWXO86ts1mo6KighUrVngTO6fTSVFR0VnPX1hYiNvt7nPB3q7XcCDOFk9XIltQUEBeXh5FRUXdkltfrkdERi4lTyIyojgcDioqKobs+Pn5+aSkpHRLLLr+Ut3bh9GysrJey3xXVlae80jKaL5Gm82GzWbr83mnlyg/V6d/eO26LpfL1SPRdLlcA6o4OJgyMzO7vbZ9JQADif/URPhM13WmY56LU6sIOp1OiouLe7xvut7bbrebwsJCcnNzWb169RnvaeuqQNhVrORUXQlVQUFBj9euPwnM2eLJzs6muLiYvLw8SktLKSkp8fl6RGRk07Q9EZFT9DY9rbS0FJvN1m0KEHQmIV0FFE7ny8jTUPP3NS5btqzHdL+u45x+/sHSNQJy+of5yspKKisrWb58+ZCc92xycnIoLS1l3bp1fS5oDAOL3+FwYLPZeoxinZ6wdiU7q1ev7nG+c11Ta/ny5RQWFnpHK/squW6z2cjOzsbpdJ61UElXIYa+KjR2LSjsyzpgfcWTk5ODy+UiPz+/138f53I9IjKyKXkSkTHp1PsmTnX6PRRut5u8vDzWrl3brV1paekZ1x7qujfj9HOeaynocxGo15ibm9trFbOBLNR7LtauXeutQldaWsqaNWtYtGgRTqdzQFO+uqoADsaH5K6kqLi4+KxV/wYS/6pVq1izZo03zq6feztmcXGx9/6e0tJScnJyBjRt9FTZ2dm43W5Wr17d4/6trpLoxcXFVFZWUlxcTGlp6VlHGwsLC8+Y3HZNmTv1fqv+6E88drsdp9NJQUFBj+T2XK9HREY4f1esEBEZLnV1dd5KYoBps9nM7OzsHhXjSkpKzJUrV5orV640MzMzu1XsqqqqMjMzM03AtNvt3ap/dZ0jOzvbtNvt3kpgXc+vqqoybTZbnxXwxtI1VlRUmCtXrjSLiorMvLy8Qa1MeKYKcRUVFabT6fRe2+nn7U+1va4Kd6dXqTub06vtdSkoKOi1yiCnVdvrT/y9XUvXe6CoqMh0OBzeanOnxuV0Oru17WpzLtUGu+I7vQJgVVVVt/dNb9UDT9dV0e70mE/HadUVV65c2SPurmN1VRvsbzynVyn05XpEZOQzTNM0hzlfExEZ09asWXPGKVqjwVi4RhkbiouLWb169ZDehygiI4em7YmIiIj0oaCgwG/3xIlI4NHIk4jIMKqsrPSuVzNajYVrlNHN7XZTW1vrXTagrq7O3yGJSIDQyJOIyDAqLy8f9UnFWLhGGd3Ky8tJSUlhxYoVfa69JSJjk0aeRERERERE+kEjTyIiIiIiIv2g5ElERERERKQfgvwdgL94PB4OHDhAZGQkhmH4OxwREREREfET0zRpbGxk8uTJWCx9jy+N2eTpwIEDJCUl+TsMEREREREJEHv37mXq1Kl97h+zyVNkZCTQ+QJFRUX5OZqxw+PxcPToURISEs6Y1cvQUR8EBvVDYFA/BAb1Q2BQPwQG9YN/NDQ0kJSU5M0R+jJmk6euqXpRUVFKnoaRx+OhubmZqKgo/ULwE/VBYFA/BAb1Q2BQPwQG9UNgUD/419lu51GPiIiIiIiI9IOSJxERERERkX5Q8iQiIiIiItIPY/aeJxEREREZHB0dHbS1tfk7jFHB4/HQ1tZGc3Oz7nkaJMHBwVit1kE5lpInERERETknpml
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1000x600 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# For the first subset\n",
|
||
|
|
"subset_indices = range(0, 13)\n",
|
||
|
|
"x_subset = magDf.iloc[subset_indices][\"GapHeight [mm]\"]\n",
|
||
|
|
"y_actual_subset = magDf.iloc[subset_indices][\"YokeForce.Force_z [newton]\"]\n",
|
||
|
|
"y_pred_subset = force_pred[subset_indices]\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot actual vs predicted for this specific condition\n",
|
||
|
|
"plt.figure(figsize=(10, 6))\n",
|
||
|
|
"plt.plot(x_subset, y_actual_subset, 'o-', label='Actual (Ansys)', linewidth=2, markersize=8)\n",
|
||
|
|
"plt.plot(x_subset, y_pred_subset, 's--', label='Predicted (Model)', linewidth=2, markersize=6)\n",
|
||
|
|
"plt.title(r\"$-15A, -15A, 0^\\circ$ roll: Model vs Ansys\", fontsize=14)\n",
|
||
|
|
"plt.xlabel(\"Gap Height (mm)\", fontsize=12)\n",
|
||
|
|
"plt.ylabel(\"Force on magnetic yoke (N)\", fontsize=12)\n",
|
||
|
|
"plt.legend(fontsize=11)\n",
|
||
|
|
"plt.grid(True, alpha=0.3)\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "b46ce772",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Better Interpolation Methods\n",
|
||
|
|
"\n",
|
||
|
|
"Let's improve interpolation accuracy without manually guessing polynomial order."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "a019f402",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Method 1: Cross-Validation to Find Best Polynomial Degree\n",
|
||
|
|
"\n",
|
||
|
|
"Automatically test different polynomial degrees and find the one with best validation score."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "fd12d56b",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Degree 1: Train R² = 0.752592, Test R² = 0.757407\n",
|
||
|
|
"Degree 2: Train R² = 0.965022, Test R² = 0.965732\n",
|
||
|
|
"Degree 3: Train R² = 0.993618, Test R² = 0.993764\n",
|
||
|
|
"Degree 4: Train R² = 0.998598, Test R² = 0.998637\n",
|
||
|
|
"Degree 5: Train R² = 0.999764, Test R² = 0.999752\n",
|
||
|
|
"Degree 6: Train R² = 0.999932, Test R² = 0.999927\n",
|
||
|
|
"Degree 7: Train R² = 0.999950, Test R² = 0.999947\n",
|
||
|
|
"Degree 8: Train R² = 0.999889, Test R² = 0.999802\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1EAAAIiCAYAAAApREJxAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAaO1JREFUeJzt/Xl0G/ed53t/iotIbWQREulVkQla8hLbkgEqixPLiQUmvfekDUp9586ZmWe6RaS388zxnZDR7eeeHs8zpxlyejlzuztzSaWX6Z47T0tAp7vvzPR0m3Baip3EjgjESZxVYsn7IlkgRNqSKJJVzx8USgBXgFsViPfrHB0DhULVF/iBND6sX33LcBzHEQAAAACgKFVeFwAAAAAA5YQQBQAAAAAlIEQBAAAAQAkIUQAAAABQAkIUAAAAAJSAEAUAAAAAJSBEAQAAAEAJCFEAAAAAUAJCFABPWJYlwzBkGIbS6fSi67a1tckwDPX09KzKvpPJpAzDkGVZRT+nv79fTU1NRa2bTqfd15b/LxwOq7+/f7llL6inp0dNTU0yDGNNto/ircXYl/LZ85vBwUF1dnaW9JxSX+9C73lTU5NisZiy2WyJVQPA0ghRADx34sSJBR9Lp9MlhR0/GRgY0OjoqEZHRzUyMqJYLKbe3l6Fw+FV20csFlMymVQ8HlcqlVIkElm1bWP51mPsy0EqlVIikViXfc1+z/v6+pRMJtXa2rrkH2oAoFSEKACeikQiGhwcXPDxEydOlG0wCAQCMk1TpmkqGAyqq6tLqVRK6XR61Y6qDQ4O6tixY4pEIgqFQgqFQquyXazMeox9ORgYGJDjOOuyr/ne85GREbW3tyscDnNECsCqIkQB8FRnZ6ey2aySyeS8jy9nOpCfBYNBRaPRBV9vKXJfCk3TXPG2sPZWc+xRvHg8LkkVFV4BrD1CFABPBYNBhUIh94tOvmQyqWw2q8OHD3tQ2drir+KVi7FfX6ZpqqurS4ODg7z3AFYNIQqA52Kx2LxT+uLxuCKRyIJHWtLptDo6OtTU1KS2trYFT9wfHBxUW1ubmpqa1NnZueA5VolEQuFw2G0EsBZHDHJH3WZPUVxs3/39/ero6JA08141NTXpt37rt9yT7zs6OmQYRsF7uNR7M982s9ms+vv7FQ6HlU6n59STzWbV2dnpbnP2mFmW5T6ee97sc1Hyt5+re771Zr+GpqYmdXR0FKy33PHKNeKYbXbDkWQy6W5/vv2XaqGxL/ZzXErty32fF/qsLPczMbtJRDGfkbWQ+6zP/tlf6jOU+91hGIZisZh6enrc9aWFf46K2fZ6/L4BsIYcAPDAyMiII8kZGhpyRkdHHUlOPB4vWEeSMzAw4N7u7u52HxsaGnKXpVIpJx6PO8Fg0IlGowXbGBgYcCQ50WjUXc80TUeSMzIyMme9gYEBJ5VKOX19fQXr9PX1OaZpFvXaUqnUnNczOjrqDA0NOcFg0DFN0xkdHS1p38Fg0AmFQk4wGHT6+voK9jMwMFCwvWLem4W22d3d7UhyQqGQMzQ0VFBzKBRy4vG4MzQ05EQikTnvYV9fn9PV1eUMDQ05IyMjTjQanfNaF9v+fO9hNBp11+vu7nY/A0u9Z8WMz9DQUMHyaDTqhEIhd7wkOX19fc7IyIi7/9nPWWjbxY59sWOVe39yPzezf1bya3ec4t/nYva/0s9E/j6L+YyU8rO20Hu+0Dq53yeOs/RnKB6Pu5+TVCrlRCIRJxgMOiMjI0v+bC617ZV8fgH4AyEKgCfyQ5TjOE4kEin44pb7ApP7cjU7RAWDwYL7+dvM/zJlmuacYJXbdu4LS+4Lc/4XLMeZ+fLY1dXlOM7yQtTsf7la8r8wFrtvSU4kEilYZ6HwWcx7s9A2c1+Y88NC7gtf/jaL+eI632tbbPv578t8gXix7ea2nXvPlhIMBuesmwtNjnMzXJSqlLHP1VHMWOV/9iKRSEFgytW+3Pd5qf2v5DOx1M/NfGO5FiEq95py+ynmMxSJRAo+I7n9zA58s3+Oltr2anx+AXiP6XwAfCEWixW0Qj5x4oSi0ei8U/lybc9jsVjB8tz5VbmW6ZZlKZvNzllvtuHhYbeG/OvM9Pf3u48tR19fn9tyOfcvHo8XvKZS9j0wMLDkPot9b5baZnt7e8FzpZtToiS5XQAzmcyCteRe58jIyJzH8qe05baf25ZlWfO+hpzVGK9oNKqTJ0+693Ofva6uLkk3X384HNbg4GDJbfaLGftSxyonFosVtP7PTaGb79zBxd7nUve/Gp+J2Rb7jKym3JTBXN3FfIZKeR35P0dLbXutft8AWF+EKAC+EI1GJd38MptIJHTkyJF518190QgEAnMeCwaD7pfL3H9zX5yWMjo6KmfmCL37L5VKlfZC8uTaLef+rXTfxbyOYt+bpbY5X73zbXO2ZDKpzs5O9xy05Sh23FYyXrmLsObOQ8m10s+9btM03W3FYjG1tbWpo6Oj6MYExYx9qWOVk/vjQu6LezweX/APDospdf/L/UzMthqfkVKdOXNGUmEQlBb/DMViMZ08edI996u3t1ehUGje92G+z+pSn8/V/n0DYH0RogD4RjQa1cDAgBukcsFqttwXlvm+ZFqW5T6e+4K31FGExba31ordd7FfkIt9bxbbZqnLczo7O9XZ2amOjg4NDQ1pdHR0WdtZ6j1ZjfGa3RUykUjMaaUfCoWUSqU0OjqqgYEBDQ8Pq7e3d9n7nK8Gqbixmq2rq8v9OUkmk/MetVvJ+7zSz8pCiv2MrKZsNqvBwcGCoFnMZ2hkZETBYNBthJHNZvXMM8/MWW/2e7Aen18A3iNEAfCNWCymZDKpEydOuNOq5pM7YjB7Klo6nVY6nXaPYOX+ajy7ffrsaTq5L9TzfUFe65bIq73vYt+b1ZbNZpVIJHT8+HF1dXUVffRvPsFgUMFgcN6phtlsdtXesyNHjujkyZNuGFmolX6uRXYkElnVTnIrGatYLCbLstTf3y/TNJd1Qer1/qys5mekFEePHlU2m1VfX5+7rJjPUDqddqdlOo6joaGhokLjUtv28vcNgNVT43UBAJCT+1KXSCQ0NDS06LrHjx93jxzk2pb39PQoEokUHME6duyYenp6ZJqmjhw5ouHh4Xkvunn8+HGFw2HFYjHFYjFlMhnF43H3v2tptfdd7HuzmnLT1np7e2WapgKBwIqO2gwMDKijo0OxWMy9IPPQ0JCGh4eVSqVW5T3r6upST0+Pent750yHSyQS6u3t1bFjx9ypbclkUseOHVv2a5rPcscqGAwqEoloYGBg0T84rNX+l2O1PyOzZTIZN4RkMhml02n19vbKsiwNDQ3NCW1LfYaCwaB6enoUi8Xco9q5ALSUpbbt5e8bAKuDI1EAfKWvr8/9q/9iotGoUqmULMtSR0eH+vr6dOzYsTnhq7u7W319fRocHNShQ4fcL+ChUKjgfI5QKKSRkRFZlqVDhw65XyyPHz+++i9yltXed7HvzWo7fvy4ex2go0ePqqOjQ5FIRG1tbSVvKxKJuK8ht738L5ir8Z7ljuDMd9QlFAqpvb3dvS5QT0+Purq61N3dXfJrWcxKxip3NGqpxilrtf/lWM3PyGy56zTlrll19OhRBYNBpVKpeX+fLPUZMk1T6XTaDfKdnZ0Kh8Nqa2tbcireUtv28vcNgNVhOI7jeF0EAAAoTe5oGc0IVl/uwsJDQ0MFASydTquzs9M9CgigcnEkCgCAMjQwMLBm57hVuuHh4XnPNQuFQopGo7QiB8CRKAAAykU2m3XP9zl69Oi6dLerRNlsVq2trYpEIorFYmpvb3fPi+vp6VnxuWgAyh+NJQAAKBPDw8Pq6OiYt+skVo9pmjp//rzbWMKyLJmmqfb29jlT/ABUJo5EAQAAAEAJOCcKAAAAAEpAiAIAAACAEhCiAAAAAKAEFd9YwrZtvfnmm9q+fbsMw/C6HAAAAAAecRxH4+Pjuv3221VVtfDxpooPUW+++aZ27drldRkAAAAAfOK1117TnXfeueDjFR+itm/fLmnmjWpoaPC0Ftu2dfHiRTU3Ny+afLF+GBN/YTz8hzHxF8bDfxg
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1000x600 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"\n",
|
||
|
|
"Best polynomial degree: 7\n",
|
||
|
|
"Best test R² score: 0.999947\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"from sklearn.model_selection import cross_val_score\n",
|
||
|
|
"from sklearn.pipeline import Pipeline\n",
|
||
|
|
"\n",
|
||
|
|
"from sklearn.model_selection import train_test_split\n",
|
||
|
|
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n",
|
||
|
|
"\n",
|
||
|
|
"# Test polynomial degrees from 1 to 8\n",
|
||
|
|
"degrees = range(1, 9)\n",
|
||
|
|
"train_scores_force = []\n",
|
||
|
|
"test_scores_force = []\n",
|
||
|
|
"train_scores_torque = []\n",
|
||
|
|
"test_scores_torque = []\n",
|
||
|
|
"\n",
|
||
|
|
"for degree in degrees:\n",
|
||
|
|
" pipe = Pipeline([\n",
|
||
|
|
" ('poly', PolynomialFeatures(degree=degree)),\n",
|
||
|
|
" ('model', LinearRegression())\n",
|
||
|
|
" ])\n",
|
||
|
|
" \n",
|
||
|
|
" pipe.fit(X_train, y_train)\n",
|
||
|
|
" \n",
|
||
|
|
" train_score = pipe.score(X_train, y_train)\n",
|
||
|
|
" test_score = pipe.score(X_test, y_test)\n",
|
||
|
|
" \n",
|
||
|
|
" train_scores_force.append(train_score)\n",
|
||
|
|
" test_scores_force.append(test_score)\n",
|
||
|
|
" \n",
|
||
|
|
" print(f\"Degree {degree}: Train R² = {train_score:.6f}, Test R² = {test_score:.6f}\")\n",
|
||
|
|
"\n",
|
||
|
|
"plt.figure(figsize=(10, 6))\n",
|
||
|
|
"plt.plot(degrees, train_scores_force, 'o-', label='Training Score', linewidth=2, markersize=8)\n",
|
||
|
|
"plt.plot(degrees, test_scores_force, 's-', label='Test Score', linewidth=2, markersize=8)\n",
|
||
|
|
"plt.xlabel('Polynomial Degree', fontsize=12)\n",
|
||
|
|
"plt.ylabel('R² Score', fontsize=12)\n",
|
||
|
|
"plt.title('Model Performance vs Polynomial Degree', fontsize=14)\n",
|
||
|
|
"plt.legend(fontsize=11)\n",
|
||
|
|
"plt.grid(True, alpha=0.3)\n",
|
||
|
|
"plt.xticks(degrees)\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"best_degree = degrees[np.argmax(test_scores_force)]\n",
|
||
|
|
"print(f\"\\nBest polynomial degree: {best_degree}\")\n",
|
||
|
|
"print(f\"Best test R² score: {max(test_scores_force):.6f}\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "86aaa7f5",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Visualize Best Polynomial Function vs Ansys Data\n",
|
||
|
|
"\n",
|
||
|
|
"Plot the continuous polynomial function to check for overfitting/oscillations between data points."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "8270efa4",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"best_degree = 6 # Degree 7 actually overfits to the data points and struggles to generalize"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "d94aa4c1",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Using polynomial degree: 6\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABKYAAAKyCAYAAADvidZRAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs/XtcXPd5L/p/1gwSAgkYQAgJSUQMkRU3dhqDnEvTNE0EubVNumOwmrP32b3Egux6J2fXJxHRrl9pfdJEgbZpdrLdFpScnN2TX1IJkpzttN2NQbk0TVrHgJ24iePIjBBC6IIEA+iGgLV+fyizDOI2MM+a+T7z/bxfL7+sGWD4rOfh4fKd9V3jeJ7ngYiIiIiIiIiIKM1CmQ5ARERERERERER24sIUERERERERERFlBBemiIiIiIiIiIgoI7gwRUREREREREREGcGFKSIiIiIiIiIiygguTBERERERERERUUZwYYqIiIiIiIiIiDKCC1NERERERERERJQRXJgiIiIiMlRXV1emIxARAQB6enoQi8UyHYOIshAXpoiI7tDf3w/HcRb9V1xcjObmZsTj8UxHNM5yNautrUVbW9uaH6+trQ3FxcUBJA1eR0cHGhsb1/xxph5zbW0tHMdBf39/pqOsi+Z5bm5uxtNPP53pGABu/0HqOA46OjoyHUWU9q/v+Vb6PtzS0mL013qyTOvXnd+3Tf0+vlb9/f2or69HbW3tgpkvKSlBfX19VnwtEZFZuDBFRLSM9vZ2jI+PY3x8HAMDA2htbUVPTw+qqqqM+aVYWn9/PxobG1FcXIzq6mo0Nzev6ePvrFlzczOOHj2K2tragBKbp6+vL2vOconFYujv70ckEkF7e3um46RE2zx3dHSgp6cHra2tmY5ijFgshra2tmW/L7W1tS25MFNdXb3s42XL1/d8d36tHzlyxP9aT+fZLiv1a629SjxeNvYrVfF4HMXFxXAcJ6n3j8ViqK+v9+u91GLzoUOH0N3djb6+PnR2dvoLUTU1NWhoaFjXky9ERCvhwhQR0TJKSkoQiUQQiUQQjUbR1NSEgYEB7N+/H7W1tVn3jGFXVxdqa2sRjUZx8uRJdHZ2or6+fk2PsVTN+vr60N/fj5aWloCSm6W9vR2e52U6hoiuri7U1dXhwQcfxIkTJzIdJyVBzXNXV1cg3wtaWlr4x/fPJc7WSpz5MzY2tuz7RiIR9PX1Lfivu7t7yffNpq/v+e78Wm9oaEBfXx+i0WhaFhSS7ddaegVkb79S1dLSgpKSkqTeNx6Po7q6GpFIBN3d3WhpaUFzc/OazmxubW1Fb28venp61huZiGgRLkwREa1RZ2cnAGTVQks8HkdjYyPa29vR2tqKmpoa/5nRVCX+MOIvsfq0t7ejsbERjY2NiMfjWdnDVOe5sbERvb29kpH8Mxjq6upEH1eruro6eJ6H8fFxRKPRVd8/8f0r8d9yH2PD1/d8dXV1aTk7cC39SrZXgH39SkZ/fz86OjqSPru5paUFNTU16OzsRF1dHZqamtDa2rro+9+xY8dQW1uL+vp6NDc3IxKJLHj7gw8+mFW/AxFR5nFhiohojSKRCJqamtDR0ZE1Z021tLT4Z5EEJVtqZYv+/n7EYjHU1dX5CySJRZxsYuI8t7a24sEHH8x0jKxmy9f3fD09PUkt7JnIxn4l49ChQzh8+PCihaPl9PT04ODBgwvuS/zcn7/QV1NT45/BttQTVI2NjX5PiIgkcGGKiGgdElvc7vylLLEdLrGN4c5ndDs6OlBdXQ3HcdDc3IyWlhb//YHb191IPHZzczOKi4v9P5ZXe+zV3r6SEydOoK6uDl1dXaivr0dxcbHYBU4Tz2zPP/sjcWHVxLWsVtpG0NLSsuTFZBPbRRI9aGtrQ21trf/YiTosdYbAap9//mPdWdPE2WWJj73z+hxLXfw2Fov5H7NSrmQkW4+enh4/e6Kfa/mcx48fRzQa9f+QraurW/bC18nWfrlM6Tqm5Sw1z6v1rLGx0Z/bxDHPv8ZLKj1PXANmKfO/du+swVJ1vLOGwPLfZ1b6/qNJIvtK31vW8vW9Hib1KXHtwP7+fuO2hybTKyDYfq1Uz7X8rEq3jo4OxGIxHDlyJKn3j8fjiMViixYnE9s+1/K9NPHz3MTr8xGRTlyYIiJah8QvdvO38CReja25uRl9fX04ePAg6uvr/T80urq60NzcjPb2dvT19SEWi6GrqwudnZ0YGBjwHycWi/mLIEeOHEEkEln1sVd7+2oSi0fHjx9HS0sLjh07ht7eXhw4cGDdNUo8ZuLC54mLOCfuq6mpwcmTJ9Ha2upv0VhK4pXT7rygeHt7+4KtH1euXEF/fz8OHTqElpYWdHd3Ix6PLzqGZD7//MdqbW31H6uxsREHDhzAwYMH0dnZiWg0iubm5lXr3NXVhZKSEr/X0WgUBw4cWNcf/QcPHlxyG8v8esTjcdTX1+PgwYMYGBhAZ2cnampqVrw2z506OjoWPFOeqM9SF3ZPpvYrZUqmxxLHtJyl5nm1nh07dgx9fX0Abp+5kbjYdLIfv5xEX2tqaha9LbFQGolE0NnZ6dfg+PHjaz7mpb7PrHS/Fok/vjs7O9HS0uL/d6e1fH2vVab7lFg0nf+qfPF4HH19fUZtD022V0Cw/QKWrudaf1alUzweR0tLy5pmNPG9cqn3LykpWfB7SDKi0eiK1wQjIloTj4iIFujr6/MAeJ2dncu+z8DAgAfAa29v9zzP88bHxxfcTjh8+LDX1NTkeZ7n1dXV+f+e/3nGx8f9+1pbWz0AXl1dnX/fao+dzOdeSeJYampqFtzf3d3tAfC6u7tXfYzEsdz5XyQS8RoaGhYcYzQa9Q4fPrxkhkTNW1tbvUgk4r+9rq5uUb47j/nw4cOL8ra3ty+qcTKff6XHmv+xS32t3Jl9KUv1LJmPm38Md/YWgNfa2up53ku9W6/EcfX19S3K3NDQsOj9k6n9aplW6/F6j2k987yUpXqWuC+ZGVluTu+UqNtSotHokvVPOHz48KKvoUTdBgYG/PuW+j6z0v3zH2e1/E1NTV5dXV3S/63Ul6WsVIP29vZF2RNfm/O/B6z163utMtWnxHG1trZ6AwMD3sDAgH/88x93vkz1K9lezT+uoPq1XD3X87Nqqe/jQdS4qanJi0aj/u2Vvm8kJOq41PerpX6mrKampkak/kREnud5OWIrXEREFkmcvn7nmRbNzc2LLkKaOPNhbGws6et7zN9usdpjJ/O5k3HndSfmn6qf7LPsra2tC65Tdeczs4lrUtyZMxqN+mcULHU9i+bmZjQ2NvrbEBJbOJa6Ds/8rIl6j42N+VsV1vL59+/fv+ix5m+xmt/btUjUZa3PUCc0NDSgo6PD/zpJnDWQqH0id21tLZqbm1FXV7ema8u0t7f7r+g1/wyfmpqaFc9QWKn2q2VarcepHtNK7pznpaTas2Q/frkzqmKxGGKxmOhWrOUeK5XPkcmtYk1NTYuuk1dfX4+2tjb09vb6X5/r/fpOhgl9mr/lrbW1FR0dHWhpaVnymkyZ6leyvQKC7dd882ux3p9Vqz2uhMQFz6XPVlrr2ZElJSW8xhQRieFWPiKidXj66acBLFy4AIDx8XF4nrfgv8R2n+bmZpw4ccK/TtHRo0dRU1Oz5C+DS/2BvNJjJ/P25SReZnq5X0rX8od44loVif/ulFhEW+qlraPR6LK/5DY0NCASifi/4Hd2dvr3rcVaP/9yWx7Wo6enB42Njaiurl7yekprkdj6ltj2dfz4cdTV1fl5Ey/Dnnjf6urqNV0z7MSJE4jH4/71cRL/JRZw1vPH4GqZVutxqse0kuXmOdWeSfY88bUpefHq5R4rExfITlw7bP5/EtevSczr/NkO4us7wcQ+tba2oqurS3QRIYh+LdUrINh+zTe/nuv9WSVpuRonvu/Pv7ZdYgEtcf29pSSOb6nvmes5Hm3bfInIbFyYIiJao3g87l/vIvGLWeIXvpV+uUtcZyZxQeR4PI6TJ08uer87f9lb7bGT+dwrSTwTvdwCVHV19boedykrZV3qoqzzNTU1+X+A9PT0LPny2Kv9ory
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1200x700 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Checking for overfitting indicators:\n",
|
||
|
|
"Number of actual data points in this condition: 13\n",
|
||
|
|
"Gap spacing: 1.50 mm average\n",
|
||
|
|
"Polynomial curve resolution: 500 points over 25.0 mm range\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was written by AI - see [5]\n",
|
||
|
|
"# Train best degree polynomial on full dataset\n",
|
||
|
|
"poly_best = PolynomialFeatures(degree=best_degree)\n",
|
||
|
|
"X_poly_best = poly_best.fit_transform(X)\n",
|
||
|
|
"model_best = LinearRegression()\n",
|
||
|
|
"model_best.fit(X_poly_best, y)\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Using polynomial degree: {best_degree}\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Select a specific operating condition to visualize\n",
|
||
|
|
"# Using the same condition as earlier: -15A, -15A, 0° roll\n",
|
||
|
|
"test_condition = magDf.iloc[0:13].copy()\n",
|
||
|
|
"currL_val = test_condition['currL [A]'].iloc[0]\n",
|
||
|
|
"currR_val = test_condition['currR [A]'].iloc[0]\n",
|
||
|
|
"roll_val = test_condition['rollDeg [deg]'].iloc[0]\n",
|
||
|
|
"\n",
|
||
|
|
"# Create very fine grid for gap heights to see polynomial behavior between points\n",
|
||
|
|
"gap_very_fine = np.linspace(5, 30, 500) # 500 points for smooth curve\n",
|
||
|
|
"X_fine = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [currL_val] * 500,\n",
|
||
|
|
" 'currR [A]': [currR_val] * 500,\n",
|
||
|
|
" 'rollDeg [deg]': [roll_val] * 500,\n",
|
||
|
|
" 'invGap': 1/gap_very_fine # Use inverse gap height for modeling\n",
|
||
|
|
"})\n",
|
||
|
|
"X_fine_poly = poly_best.transform(X_fine)\n",
|
||
|
|
"y_fine_pred = model_best.predict(X_fine_poly)\n",
|
||
|
|
"\n",
|
||
|
|
"# Extract actual Ansys data for this condition\n",
|
||
|
|
"gap_actual = test_condition['GapHeight [mm]'].values\n",
|
||
|
|
"force_actual = test_condition['YokeForce.Force_z [newton]'].values\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot\n",
|
||
|
|
"fig, ax = plt.subplots(1, 1, figsize=(12, 7))\n",
|
||
|
|
"\n",
|
||
|
|
"ax.plot(gap_very_fine, y_fine_pred[:, 0], '-', linewidth=2.5, \n",
|
||
|
|
" label=f'Degree {best_degree} Polynomial Function', color='#2ca02c', alpha=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot actual Ansys data points\n",
|
||
|
|
"ax.scatter(gap_actual, force_actual, s=120, marker='o', \n",
|
||
|
|
" color='#d62728', edgecolors='black', linewidths=1.5,\n",
|
||
|
|
" label='Ansys Simulation Data', zorder=5)\n",
|
||
|
|
"\n",
|
||
|
|
"ax.set_ylabel('Force (N)', fontsize=13)\n",
|
||
|
|
"ax.set_title(f'Degree {best_degree} Polynomial vs Ansys Data (currL={currL_val}A, currR={currR_val}A, roll={roll_val}°)', \n",
|
||
|
|
" fontsize=14, fontweight='bold')\n",
|
||
|
|
"ax.legend(fontsize=12, loc='best')\n",
|
||
|
|
"ax.grid(True, alpha=0.3)\n",
|
||
|
|
"\n",
|
||
|
|
"# Add vertical lines at data points to highlight gaps\n",
|
||
|
|
"for gap_point in gap_actual:\n",
|
||
|
|
" ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.3, linewidth=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"# Check for excessive oscillations by looking at derivative\n",
|
||
|
|
"print(\"Checking for overfitting indicators:\")\n",
|
||
|
|
"print(f\"Number of actual data points in this condition: {len(gap_actual)}\")\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Gap spacing: {np.diff(gap_actual).mean():.2f} mm average\")\n",
|
||
|
|
"print(f\"Polynomial curve resolution: 500 points over {gap_very_fine.max() - gap_very_fine.min():.1f} mm range\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "8e4487e0",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Compare Multiple Conditions Side-by-Side\n",
|
||
|
|
"\n",
|
||
|
|
"Let's check several different operating conditions to ensure consistent behavior."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "74d6e98f",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAS1CAYAAADjrFn4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs/X1wW9l5J/h/L0iKeqFIEBT10t1qt8Bud2zHjg1SbsfxxC8i7XGcTe+4QfVuTXbjikfkuuza3XI2ZDSp2q3e3+wqZCZVuzX2xKR6Zu2ajackMjPOzqQzMdFx3uzYLRFuO7Zju5tQv+qFEkmQFCVRJHF/f8C4AkgCwn0u8ODg8Pup6mqBwME557kPL87hwbnXcV3XBRERERERERERERERkeFCtW4AERERERERERERERFRObioQUREREREREREREREdYGLGkREREREREREREREVBe4qEFERERERERERERERHWBixpERERERERERERERFQXuKhBRERERERERERERER1gYsaRERERERERERERERUF7ioQUREREREREREREREdYGLGkRERGS9dDoNx3G2/a+7uxvj4+O1bmJdKBbH9vZ29PX1IZlMit97fHwcjuMglUpVsMX6xsfH0d3dLYpFpWOQOz7Dw8P3fW17ezscx0F/f7+oLr9tNylOm6XTaQwODnox6erqwuDgINLpdFXq0xYk9tWSa1MuZ7u6utDf349EIqHelu7ubnR3d29pW7nxMjG+RERERLbhogYRERHtGL29vZiZmfH+m5iYQE9PDwYHB9HX11fr5qnI/cG2q6vL++Oh3z8c5sdxenoaZ8+eBZD9Y+Do6Gg1ml03ZmZmkEwmMT8/X+umeCYnJ0s+n0gk1P9gb2KcACCZTOLYsWNIJBIYGBjAxMQE4vE4zp8/b80fqk2KfTqdRl9fHwYHB9HT04OJiQmMjY0hHo8jkUhgYmKi1k30HS+T4ktERERkq8ZaN4CIiIhISzQaRTQaLXgcj8fR39+Pvr4+dHd3Y3p6uoYtrK5UKoW+vj6Ew2GMjIwgGo0ilUohEon4ep/NcYzFYojH4xgcHMTw8DB6e3sRi8Uq3fy6MDIygpGRkVo3wxOLxZBMJpFMJosek4mJCe91lTY5OYloNLqlbtPiBGQXNLq7u9Hb24uJiQmEw2EAQDwex+nTp3HixAl0d3djZmamIP9NVQ+x7+7uRiqVwtTUFHp7ewueGxkZMWLnVrF41UN8iYiIiGzFnRpERES04/X29mJsbAzJZPK+32qvZ319fYhGo5ienkY8HvcWIyq1ADE2Nlbwf6q9np4eRKPRksdkfHwcTz/9dFXqP3XqFM6dO1eV9660U6dOIRwOY2pqylvQyAmHw3j++ecBAIODgzVonX+mx354eLjogkaOyYtHpseXiIiIyGZc1CAiIiICMDAwgHA4XNb9B+rR+Pg4UqlU1S/nEg6Hjfh2NWWFw2Hv8knbyS3iDQwMaDbLOIlEAslksuQ37HM7nHKvJbl0Oo3R0VFEo9GiCxpERERERMVwUYOIiIjoZ06ePIlUKrXl/gLpdBr9/f1ob29HV1fXtgsfuUvX5N98PFcmd7+OyclJtLe3e//u7u4u+NZ3OfWU85rtjI2NIRaLYX5+Hv39/ejq6kJXV1dFb5KeTqeRTqe3fLs6/z4e7e3t6O/vv+/Cx/DwMBzH2fZeD/k38s3FdHMd232bvpx2bH6/9vb2gvfLP85dXV1bdvZsdxPrdDqN4eHhgvuY+Ln3SJAbeAPA008/jXQ6ve29U86dO4dYLLZlZ0JOsZtyb76Z8mb9/f3e8RsdHfVuAJ3Lt+3eN5FIeJcjysUrd5zKIf3dAO7tLrrfH9hzz+d/Q99vu+/XzlLniXJySRJ7P79H5Zzr7ie3yOZ3EVnyO1yqL7ljl+vLdjv1NsdLEt9Ktz13KcFc3byfEREREe00XNQgIiIi+pmuri4AwMWLF72fpVIpHDt2DKlUCmfPnsXw8DDGx8cL/hiXTqfR3d2Nnp4eTE9Po7e3F4ODgzh+/Dief/5579vf8/Pz3h/C+vv7EYlEvD8CllNPOa8pJrdY09fXh76+Pu+eGoODg75vFF5Mrh2bF2pyN14eGRnB2bNnvXiV+rZ77j3OnDmzpR/JZNJ7PhfT7u5u75v0vb29GB8fL/gjX7nt2Px+Z8+e9d6vr68P/f39OH36tPdH8P7+/vveZPv8+fNIJpMYHh7G1NSUd++Rcv+gG4/Hcfz48bJeu51YLLbtJajS6TQmJyercjmlkZER7/408Xgc09PTmJ6exsmTJ4uWyR3b3GLJ8PAwenp6MDk5ed+FjSC/G7nywP0vd5R7Pj9n/LS7nHaWOk+Uk0uS2Pv5PSrnXHc/MzMzBfEsh/R3uFhfkskk+vr6kEqlMDY25v1e328XjiS+lW57LtempqYwMTGB3t5eTE1NlR1LIiIiorrnEhEREVluYWHBBeAODAyUfN3ExIQLwB0bG/N+1tvb64bD4YLXTU1NufnDqKGhoS2viUajW+obGxtzAbgA3KmpqYLnyqmnnNdsJ9d/AO7MzMyWdsZisZLlN79Pb2+vOz097f03MTHh9vb2ugDckZGRgjLxeNyNRqNb3isWixXUm4tNfvu26+/Q0JALwF1YWCgotznWuXZK27H5/cLhsAvAnZ6e9n6Wi31+vmzXj+3E4/Etx63csuUC4A4NDbmuey9u29WXiyUANx6Pl9Wmco6f62bjlmvD/d4397OJiYmC1+ZyK/facnOlnN+NnGg0uqV8MQC27fv92l1uO0udJ7azXS65riz29/s9Kvdcdz8DAwO+cz3o7/DmvuSOTy7/czafEyuR25Vs+/T0dEXPE0RERET1iDs1iIiIiH5mfn4ewL1vD+cu2TMwMOBdWimdTqOnpwfhcNjb4ZBMJtHT01PwXtFotOglloaGhgouc1NOPeW2pZTe3t4t34yOx+NIJpP33W2QL3fJltx/p06dApD91vDQ0FBBvyYnJ7f9tvzIyAiSyeR9d2tsvmzS5OQk4vH4lsslbf5GfO5SW9J2bH6/aDSKcDhccFP1/DzxK7fzQlJWIrcbI//yOrlveBe79FStbL5xfe7b/9tdGggo//e0lHA47JUrJfc7HYlEfLfbbzs3nyeKqWQulfo9Avyf64rJ5Vy55SrxO7z5nJBIJDA0NLTtTeErqdJtz7VvZGRE7fxBREREZBouahARERH9TO6SIrk/Vuf+4DY6OurdWyH3XzqdLrhkzeY/zqVSqaKXVtl83fly6im3LdvJ/RFsu/bkLrnl54+S8Xgcrut6/y0sLGBqamrLH2Bzl/Ha/EfQ/J/lX+pru3rC4bB32aRkMolUKrXt5ZJKXcZG0o7t3m/zH7K3+8N2Mbl7IXR3d6O9vV39hvTRaBSxWMyLZe6PutW49FSl5RYLcpcs2izI70bO5t/5Ysq9TNV27fbbzmL3p6hmLpVz+S0/57pi/J53KvU7nJOrN9eOaqp026PRKAYGBjA+Pi66fwwRERGRDbioQURERPQziUQC4XB4yx+UpqenC/6In/tvYGAAQPZb8KlUquBGscX++A4U/4PV/eop9zXbCYfDBd+41pD7FvF29eZ/Y7qUgYEB75vu586dQzgcLuvb65VuByD/Bnfu5r/JZBKnT5/G9PR0wY4WLU8//bS36yeXq/F4XL0d1SL93QCysQGw5b4jm+XfS6Xa7dzuPFHrXPJ7rismd/+J+8U7p1K/w5v5WZiUqkbbx8bGMDU15eXM6Oho2TdpJyIiIrIBFzWIiIiIkP1jYSqVKrjRbe6PiqV2E+Sej0ajGB4ehuM4GB4exsTExJbL0RRTTj3ltqWYnp6ebS/1lPsWeblt9SP3ntvVm/vZ/eo9ffo0gOzxmZycLHkz3mq2I4j+/n7E43FMTU0hHo8jGo2io6OjavUVk/sDaCKRwLlz5wItaGhe9iZ3jHI3R94s6O8GkF3cicViGB8fL9q3VCqFyclJ9Pb2lrWwtrndlWhnrXMp6LkuJxwOY2hoCMlksqzLg1X6dzh3LC5cuFB2GalqnX96e3sxNjaGmZkZjIyMIJFI+L4
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1600x1200 with 4 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Visual inspection guide:\n",
|
||
|
|
"✓ Look for smooth curves between data points (good)\n",
|
||
|
|
"✗ Look for wild oscillations or wiggles between points (overfitting)\n",
|
||
|
|
"✓ Polynomial should interpolate naturally without excessive curvature\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was written by AI - see [6]\n",
|
||
|
|
"# Select 4 different conditions to visualize\n",
|
||
|
|
"# Each condition has 13 gap height measurements (one complete parameter set)\n",
|
||
|
|
"condition_indices = [\n",
|
||
|
|
" (0, 13, '-15A, -15A, 0°'), # First condition\n",
|
||
|
|
" (13, 26, '-15A, -15A, 0.5°'), # Same currents, different roll\n",
|
||
|
|
" (91, 104, '-15A, -10A, 0°'), # Different right current\n",
|
||
|
|
" (182, 195, '-15A, -10A, 1°') # Different right current and roll\n",
|
||
|
|
"]\n",
|
||
|
|
"\n",
|
||
|
|
"fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n",
|
||
|
|
"axes = axes.flatten()\n",
|
||
|
|
"\n",
|
||
|
|
"for idx, (start, end, label) in enumerate(condition_indices):\n",
|
||
|
|
" ax = axes[idx]\n",
|
||
|
|
" \n",
|
||
|
|
" # Get condition data\n",
|
||
|
|
" condition_data = magDf.iloc[start:end]\n",
|
||
|
|
" currL = condition_data['currL [A]'].iloc[0]\n",
|
||
|
|
" currR = condition_data['currR [A]'].iloc[0]\n",
|
||
|
|
" roll = condition_data['rollDeg [deg]'].iloc[0]\n",
|
||
|
|
" \n",
|
||
|
|
" # Create fine grid\n",
|
||
|
|
" gap_fine = np.linspace(5, 30, 500)\n",
|
||
|
|
" X_condition_fine = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [currL] * 500,\n",
|
||
|
|
" 'currR [A]': [currR] * 500,\n",
|
||
|
|
" 'rollDeg [deg]': [roll] * 500,\n",
|
||
|
|
" 'invGap': 1/gap_fine # Use inverse gap height for modeling\n",
|
||
|
|
" })\n",
|
||
|
|
" \n",
|
||
|
|
" # Predict with best degree polynomial\n",
|
||
|
|
" X_condition_poly = poly_best.transform(X_condition_fine)\n",
|
||
|
|
" y_condition_pred = model_best.predict(X_condition_poly)\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot polynomial curve\n",
|
||
|
|
" ax.plot(gap_fine, y_condition_pred[:, 0], '-', linewidth=2.5, \n",
|
||
|
|
" color='#2ca02c', alpha=0.8, label=f'Degree {best_degree} Polynomial')\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot actual data\n",
|
||
|
|
" ax.scatter(condition_data['GapHeight [mm]'], \n",
|
||
|
|
" condition_data['YokeForce.Force_z [newton]'],\n",
|
||
|
|
" s=100, marker='o', color='#d62728', \n",
|
||
|
|
" edgecolors='black', linewidths=1.5,\n",
|
||
|
|
" label='Ansys Data', zorder=5)\n",
|
||
|
|
" \n",
|
||
|
|
" # Formatting\n",
|
||
|
|
" ax.set_xlabel('Gap Height (mm)', fontsize=11)\n",
|
||
|
|
" ax.set_ylabel('Force (N)', fontsize=11)\n",
|
||
|
|
" ax.set_title(f'currL={currL:.0f}A, currR={currR:.0f}A, roll={roll:.0f}°', \n",
|
||
|
|
" fontsize=12, fontweight='bold')\n",
|
||
|
|
" ax.legend(fontsize=10, loc='best')\n",
|
||
|
|
" ax.grid(True, alpha=0.3)\n",
|
||
|
|
" \n",
|
||
|
|
" # Add vertical lines at data points\n",
|
||
|
|
" for gap_point in condition_data['GapHeight [mm]'].values:\n",
|
||
|
|
" ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.2, linewidth=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.suptitle(f'Degree {best_degree} Polynomial: Multiple Operating Conditions', \n",
|
||
|
|
" fontsize=15, fontweight='bold', y=1.00)\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"Visual inspection guide:\")\n",
|
||
|
|
"print(\"✓ Look for smooth curves between data points (good)\")\n",
|
||
|
|
"print(\"✗ Look for wild oscillations or wiggles between points (overfitting)\")\n",
|
||
|
|
"print(\"✓ Polynomial should interpolate naturally without excessive curvature\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "81a46ab5",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Find Worst-Performing Segment\n",
|
||
|
|
"\n",
|
||
|
|
"Let's identify which parameter set has the worst R² score to showcase model performance even in challenging cases."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "35f4438e",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Total segments analyzed: 637\n",
|
||
|
|
"\n",
|
||
|
|
"Worst performing segment:\n",
|
||
|
|
" Rows: 7631 to 7644\n",
|
||
|
|
" Parameters: currL=15A, currR=15A, roll=3.0°\n",
|
||
|
|
" R² Score: 0.999302\n",
|
||
|
|
"\n",
|
||
|
|
"Best performing segment R²: 0.999992\n",
|
||
|
|
"Mean R² across all segments: 0.999934\n",
|
||
|
|
"Median R² across all segments: 0.999953\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was written by AI - see [7]\n",
|
||
|
|
"# Calculate R² score for each 13-row segment (each parameter set)\n",
|
||
|
|
"num_segments = len(magDf) // 13\n",
|
||
|
|
"segment_scores = []\n",
|
||
|
|
"segment_info = []\n",
|
||
|
|
"\n",
|
||
|
|
"for i in range(num_segments):\n",
|
||
|
|
" start_idx = i * 13\n",
|
||
|
|
" end_idx = start_idx + 13\n",
|
||
|
|
" \n",
|
||
|
|
" # Get actual and predicted values for this segment\n",
|
||
|
|
" segment_actual = y.iloc[start_idx:end_idx, 0].values # Force only\n",
|
||
|
|
" segment_pred = model_best.predict(poly_best.transform(X.iloc[start_idx:end_idx]))[:, 0]\n",
|
||
|
|
" \n",
|
||
|
|
" # Calculate R² for this segment\n",
|
||
|
|
" segment_r2 = r2_score(segment_actual, segment_pred)\n",
|
||
|
|
" segment_scores.append(segment_r2)\n",
|
||
|
|
" \n",
|
||
|
|
" # Store segment info\n",
|
||
|
|
" segment_data = magDf.iloc[start_idx:end_idx]\n",
|
||
|
|
" currL = segment_data['currL [A]'].iloc[0]\n",
|
||
|
|
" currR = segment_data['currR [A]'].iloc[0]\n",
|
||
|
|
" roll = segment_data['rollDeg [deg]'].iloc[0]\n",
|
||
|
|
" segment_info.append({\n",
|
||
|
|
" 'start': start_idx,\n",
|
||
|
|
" 'end': end_idx,\n",
|
||
|
|
" 'currL': currL,\n",
|
||
|
|
" 'currR': currR,\n",
|
||
|
|
" 'roll': roll,\n",
|
||
|
|
" 'r2': segment_r2\n",
|
||
|
|
" })\n",
|
||
|
|
"\n",
|
||
|
|
"# Find worst performing segment\n",
|
||
|
|
"worst_idx = np.argmin(segment_scores)\n",
|
||
|
|
"worst_segment = segment_info[worst_idx]\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Total segments analyzed: {num_segments}\")\n",
|
||
|
|
"print(f\"\\nWorst performing segment:\")\n",
|
||
|
|
"print(f\" Rows: {worst_segment['start']} to {worst_segment['end']}\")\n",
|
||
|
|
"print(f\" Parameters: currL={worst_segment['currL']:.0f}A, currR={worst_segment['currR']:.0f}A, roll={worst_segment['roll']:.1f}°\")\n",
|
||
|
|
"print(f\" R² Score: {worst_segment['r2']:.6f}\")\n",
|
||
|
|
"print(f\"\\nBest performing segment R²: {max(segment_scores):.6f}\")\n",
|
||
|
|
"print(f\"Mean R² across all segments: {np.mean(segment_scores):.6f}\")\n",
|
||
|
|
"print(f\"Median R² across all segments: {np.median(segment_scores):.6f}\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "080f954b",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Visualize Worst-Performing Segment\n",
|
||
|
|
"\n",
|
||
|
|
"Even our worst case shows excellent model performance!"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "6ceaf2b2",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAKyCAYAAAAEvm1SAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3X9Y3Nd9J/r3d0BCSAIGEAghiUiDFeVHk9iDXDfpttvEg9Pbu5s2NVjNPnu3270WpPWT3NZbQ9QqbdzWUSDZPL3OehOQ02x72ybS4HTr/tjGjOM2bZI6Etipmx+KzQhjGQQIGEACITHf7/1jPF+BQDAw5zCfw3m/nsePNcP8eJ/z4TPAme/5juN5ngciIiIiIiIiIqINFMh1ACIiIiIiIiIisg8XpYiIiIiIiIiIaMNxUYqIiIiIiIiIiDYcF6WIiIiIiIiIiGjDcVGKiIiIiIiIiIg2HBeliIiIiIiIiIhow3FRioiIiIiIiIiINhwXpYiIiIiIiIiIaMNxUYqIiMhCXV1d675vPB5HfX09SktLUVpaisbGRsTjcYXpSJVYLMbaEBERkVhclCIi4zU3N8NxnNt+3XEc1NbWLvu1eDwOx3EQi8V0xVMiHo+jsbERpaWlcBwHpaWlqK+vR2dnZ66jAQB6e3vhOM6S/+rq6tDe3q78+VpbW/250PH4m11zczPOnj3rX16pfq2trUgkEovuX1tbi3A4jGg0ira2NvT29qKurm7J7Sj3ysrKUF9fv6lrE4vF4DiOmNfDtaqrq4PjOOjt7c11lKyt9bVEKmk1aW9vR2lp6arXSdXb24v6+nrU1dUZ26dEpA8XpYjIePX19QCw7MJS+rp4PL7sL8Ppr0ciEX0Bs9TV1YXa2lrE43GcOnUKPT09aGtrAwBEo9Ecp1uso6MDExMTmJiYQF9fH5qbm3Hy5EnU1dUpe47m5mbEYjFEo1H09PSIrp1EnZ2diMVi/vfQQrfW7/jx44jFYjh48OCio23Si1GRSARNTU3o7u5GIpHI6eJu+uit9CJ0pn/4LFzwra2tXfYIskxuk+ntYrGY/wdvbW3tbRdVMxlPZ2cnamtr/T/6l7tNOBxGQ0MDGhsbV5sKykA8Hkd7ezuam5uX/Xp7e/uyizIrvTHS29uLYDCIjo4OndE3VKavJaqsVBfWJDOZvjYtZ7XXq2PHjqG7uxs9PT2IRqPGLE4S0QbxiIgMNzEx4QHwWlpalnytqanJa2ho8AB40Wh0ydcbGhq8SCSiPWM0GvUmJibWfL+enh4PgNfQ0LDs19fzmDqkcy43x319fbetz3rc7nkoM8Fg0Ovu7l503Ur18zzPC4fDXjgcvu1jdnd3ewC8np4epVkzlX4NaGho8Lq7u72Ojg4PgNfW1rbi/fr6+rxgMOg1NTV5PT09/v06OjrWdJtMbxeNRr1gMOh1dHR4fX19/m2amprWPJ6Ojg4vGAx60WjU6+7u9pqamlbss+Xqvlmkv/9urYmO5wgGgyu+Jre1tXnBYNDr6elZ9F9fX99tbx+JRLympiYvGAxqy79Rsn0tWatM6rIZapIew2rXrVemr03LyeT1amHNI5GImN9diEgGLkoR0aYQDoe9UCi05PpQKOR1dHR44XB42V9Wg8Hgqn+4qgBgXX8QhkKhZcclzWp/iDQ0NCj5QyT9y+9m/eNat/RCxq1Wq19LS4t3u/exuru7vXA4nNOFwqampiXfX21tbbfNnBaJRJa938I5yuQ2a3msW793l5vbTMaT/iN7ucda7o/t5R5zs9iIRamFQqHQqotSa3msjo4Ofwymv7Zl81qSrdvVZTPURPeiVKavTcvJ5PWqp6fHC4fDXiQS4ZtKRLQEt+8R0aYQiUSWbNFLJBKIx+OIRCKIRCJLthb19vYikUiI3f6VPkFxa2trrqMowcP1c6+trQ0PPPDAmu8Xi8UQCoUWXZc+R0g8HkdPTw8aGhpUxVyzWCyGo0ePLrquqanJ/9py0tsNb93u09DQ4H8tk9tk+lgA0N3dndHrTSbjiUQiCIfDi26Tfv7lzoPT2NiI3t5envRckHQ90j+jAHlbslVb7rVEEhtrAmT+2rScTF6vwuEwenp60N3dndOfFUQkExeliGhTWO68UrFYDMFgEKFQyD/R78I/yNJfX/iHXfoP7fQ5YW49p0J7e7v/XM3NzSgtLfX/6EyfiyF9EvLe3l40Njb6J2FPn29hpZOyL5T+wzLTRYRbT4ZeV1e35I/T2+VM6+rq8r9eV1en5BxB6fm59RfelZ5ruXn+3d/9Xf+krum5XHjeivXWrr293Z+rW/MkEolF5wi69TwZmcz5wsdP517udreOQUd90uf9yFT6e7i3t3fRuVWam5vR2NiI5uZmhEIhxGIxf742Wrqvb/1DNxgMIhgM3vZExePj4wBSJwJfKP04vb29Gd0m08e6na6urkWvQZmOZ6U/lNN5Fkr3n+4TN6/0PZz+gIKF0icpX/javFKvLne9ROl8K52b5/Tp0wiFQn6tI5GI0pNAS6rF7V5LNlKua7LSnK328ysXbn1tWs56X3+JiBbiohQRbQrpP7i6u7v9606fPu0v6KS/vvCP+FvfGUwv2ITDYTz77LNoa2tDR0fHkhMEx+Nxf0Hg+PHjAFKLJEePHkVfXx+i0SjC4TDGx8f9E5MDqT8i0yd+zUT609GCwWBGt+/q6kJZWRmi0Sj6+voQCoVw7733+r/0JhKJ2+YEUidNTi809PT04OjRo/6RMOuxcLEOwKITa2fyXLfO86OPPurPZfokugvfjV1P7YLBIMbGxtDb24tjx46hra3NP2l3Y2Mj7r33Xhw9ehTRaBShUAjNzc2LMq425wAWPX5ra6v/+Pfee++ibOlFsWAwiGg06tfn9OnTSuqz8B3r20kvoi78xKxEIrHohPKJRAKdnZ3+glx9fb3/35kzZzLKolL6+3e5PikrK0NfX9+y90svIN26gJOu3djYWEa3yfSxForH44t649lnn816PMDqH9wQCoUWvUaqttr38Fos16srXS9J+g/1aDSK1tZW/79bdXZ2LjpqJP16dbsT6a9FrmuRyWvJRpJQE2D5Ocv059dGWOm1aTnZvF4REaXl5zoAEZEqt27Ri8ViOHXqlH85HA6ju7t70ULGrUd/tLS0+Isn4XAY4XDY/xSt9C+q6UP70wsk6edsaWkBkPrDb+Ev3el3ENPvHGZqrUcApJ8/7dSpUygtLcWZM2fQ1NSEc+fO3TZnIpFAc3MzOjo6/PkJh8MYGxvzf0HOxK2/RAeDQUQiEZw6dcofe6bPdes8pzMDqV92F87lemu3UPrT5IDUUQTNzc2IRCL+fcvKyvwjnNI5Vpvz1R4/kUj442hsbERDQ8OiI2BU1ie9eLXS1pm2tjZ/vB0dHWhvb0dHR8ei+wSDQXiet+rz3c6tC3uZ3H6l7R7rPVImfRTlwtcEAP7CWro2q90m08daaOFiYkdHx6Lv5WyO/Glra0NLS8ttaxwMBpc9ikqVlb6H1+p2vbpSD69G9ffectKveQsX//r6+tDe3r5o4Sa9fXzhtqcHHngAzc3NOH36dNZbnHJdi0xeS9J010VKTYDl5yzTn1+rUTGPK702LUfqkYpEZBYuShHRplFfX4/W1lYkEgmMj48vOV/UwsPwb90alz6PxK3nhAmFQv67ywt/eVu4CHDkyBEAQF1dnb+QoeKcGenHWO7Q+Eykf5lMv1O5Us70glVzc/OSOVjt8P2F2traFv1RvtwvtGt5rkwWW7Kp3ULp+UnfF8CirW7pbCv9UX/rnC+03ELl+Pg4gsEg4vE44vH4bbOpqE8mfzws3LbS1taGzs5OtLa2Kj2nykZv3Vnpj6r04mBrayuOHj2Kc+fO+fnSHxefyW3Wcjvg5vdHektTNBrN+Aim242nvr4ekUhk0RGJtyorK9N2TqnVvofX43aPtd7n2IjvvaampiUL0vX19Whvb8e5c+f814H0H/yhUGhRb4bD4ayPypFQi7W8luiui4SaLLRwvGv9+ZXp465XNq9Ny5F4JCMRycPte0S0aaR/cUuf3yYUCi36hejo0aNIJBLo7e1d8vX
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1200x700 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Even in the worst case, the model achieves R² = 0.999302\n",
|
||
|
|
"This demonstrates excellent interpolation performance across all operating conditions!\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was written by AI - see [8]\n",
|
||
|
|
"# Get worst segment data\n",
|
||
|
|
"worst_data = magDf.iloc[worst_segment['start']:worst_segment['end']]\n",
|
||
|
|
"worst_currL = worst_segment['currL']\n",
|
||
|
|
"worst_currR = worst_segment['currR']\n",
|
||
|
|
"worst_roll = worst_segment['roll']\n",
|
||
|
|
"worst_r2 = worst_segment['r2']\n",
|
||
|
|
"\n",
|
||
|
|
"# Create fine grid for this worst segment\n",
|
||
|
|
"gap_fine_worst = np.linspace(5, 30, 500)\n",
|
||
|
|
"X_worst_fine = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [worst_currL] * 500,\n",
|
||
|
|
" 'currR [A]': [worst_currR] * 500,\n",
|
||
|
|
" 'rollDeg [deg]': [worst_roll] * 500,\n",
|
||
|
|
" 'invGap': 1/gap_fine_worst\n",
|
||
|
|
"})\n",
|
||
|
|
"\n",
|
||
|
|
"# Get predictions\n",
|
||
|
|
"X_worst_poly = poly_best.transform(X_worst_fine)\n",
|
||
|
|
"y_worst_pred = model_best.predict(X_worst_poly)\n",
|
||
|
|
"\n",
|
||
|
|
"# Extract actual data\n",
|
||
|
|
"gap_worst_actual = worst_data['GapHeight [mm]'].values\n",
|
||
|
|
"force_worst_actual = worst_data['YokeForce.Force_z [newton]'].values\n",
|
||
|
|
"\n",
|
||
|
|
"# Create the plot\n",
|
||
|
|
"fig, ax = plt.subplots(1, 1, figsize=(12, 7))\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot polynomial curve\n",
|
||
|
|
"ax.plot(gap_fine_worst, y_worst_pred[:, 0], '-', linewidth=2.5, \n",
|
||
|
|
" label=f'Degree {best_degree} Polynomial', color='#2ca02c', alpha=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot actual data points\n",
|
||
|
|
"ax.scatter(gap_worst_actual, force_worst_actual, s=120, marker='o', \n",
|
||
|
|
" color='#d62728', edgecolors='black', linewidths=1.5,\n",
|
||
|
|
" label='Ansys Data', zorder=5)\n",
|
||
|
|
"\n",
|
||
|
|
"# Formatting\n",
|
||
|
|
"ax.set_xlabel('Gap Height (mm)', fontsize=13)\n",
|
||
|
|
"ax.set_ylabel('Force (N)', fontsize=13)\n",
|
||
|
|
"ax.set_title(f'Worst Case Performance (R² = {worst_r2:.6f}) - currL={worst_currL:.0f}A, currR={worst_currR:.0f}A, roll={worst_roll:.1f}°', \n",
|
||
|
|
" fontsize=14, fontweight='bold')\n",
|
||
|
|
"ax.legend(fontsize=12, loc='best')\n",
|
||
|
|
"ax.grid(True, alpha=0.3)\n",
|
||
|
|
"\n",
|
||
|
|
"# Add vertical lines at data points\n",
|
||
|
|
"for gap_point in gap_worst_actual:\n",
|
||
|
|
" ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.3, linewidth=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Even in the worst case, the model achieves R² = {worst_r2:.6f}\")\n",
|
||
|
|
"print(f\"This demonstrates excellent interpolation performance across all operating conditions!\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "a8f24086",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Torque Prediction Analysis\n",
|
||
|
|
"\n",
|
||
|
|
"Now let's examine torque predictions. Torque is influenced by multiple variables:\n",
|
||
|
|
"- **Roll angle** (primary driver of torque asymmetry)\n",
|
||
|
|
"- **Current imbalance** (currL vs currR difference)\n",
|
||
|
|
"- **Gap height** (via invGap feature)\n",
|
||
|
|
"\n",
|
||
|
|
"We'll visualize torque accuracy by sweeping across roll angles first, then examining current variations."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "d896f5a7",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Overall Torque Prediction Accuracy\n",
|
||
|
|
"\n",
|
||
|
|
"First, let's check the overall R² score for torque predictions across all data."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "c120048d",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAMWCAYAAADs4eXxAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAuftJREFUeJzs3XlcVPX+x/H3DLsIDKioqamjppVWjli3vXRsXxX0ttzq3ptgd2kP1LpLdUux/XZbwOr+bt26KbTvibbZKqCV2aKMliI6JAwzKoPCzO8PmpEBVHSAgeH1fDx4PJrzPef4mXMG4s33fL9fg9fr9QoAAAAAALQ7Y6gLAAAAAAAgXBG6AQAAAADoIIRuAAAAAAA6CKEbAAAAAIAOQugGAAAAAKCDELoBAAAAAOgghG4AAAAAADoIoRsAAAAAgA5C6AYAAGiisLDwoI+12WyaMmWKkpOTlZycrIyMDNlstnasDu2lqKiIewOgUxC6AaATFRUV+X8Zb8tXT/iFsLS0VAaDocXXhAkTtGDBgg79t4uKimQwGAKuc35+vjIyMjr0392f1urqCiZMmCCDwaDS0tJQl7JfCxYsUHJy8gEfl5WVpRUrVvhf7+vzmZOTI4fDEXD8iBEjZLFYVFBQoNzcXJWWlmrChAkt9kPopaSkaMqUKdwbAB2O0A0AnSgtLU0FBQUBX9OnT5fD4dCcOXNatJnN5lCX3Gny8vJUXV2t6upqlZWVKSsrS/PmzdOECRM6tY6SkpKgejrDlc1mU2lpqUwmk/Ly8kJdTofIz89XUVGRcnNzW7Q1/3zOmTNHRUVFGj58eMAfR3xh22q1KjMzU0uWLJHD4VBRUVFnvpUAvt53g8GgESNGKD8/v83HZWRkKDk5WSNGjGj1+6It+7R1v9LS0oA69/ZHt/2da8GCBa3+oWTEiBEB+1ksFqWnp4f8j2wAwh+hGwA6kclkktVqDfjyhcrm261Wa4ir7VwpKSkymUwymUwym83KzMxUSUmJSktLlZOT02l15OXlyev1HvBxhYWFYd1jVlhYKKvVqunTp2vx4sVBn6srXqucnJy9/kGh+eczPT1dJSUlMpvNAaEtPT094DhfIA/VH9AcDodGjBghk8mkJUuWKCcnR1lZWft9isRms2nChAlKSUnR0qVLlZOTo4yMjIDA3pZ92rqf74mAKVOmqKysTHl5ecrLy2sRiNv6b5pMJpWUlAR8LVmypMX7zM3NVXFxcUj/KAKgB/ACAEIqLy/PK8lbUlIS6lJCoqSkxCvJW1BQ0Gp7enq612KxdMi/vWTJEq8kb1lZWdDnkuRdsmRJO1TVvnW1F7PZ7M3Ly/PXFsx7bc9rtTe5ublek8nU5v3z8vJa3X9/n8/s7Gzv3n6dWrJkiddisez12M6QmZnZ4vsnNzd3rzX7WK3WVo9reo3ask9b97NYLN7MzMyAfXzXvun3QVvOdaD3vrVrBADtiZ5uAECX1xV7RXuS0tJS2Wy2gCcwCgoKQlxV+8rNzdX06dMP+LiioqIWvdi+x6RtNptKSkpa9H53pqKiIs2YMSNgW2Zmpr+tNb7H4bOysgK2p6en+9vask9bzyXt6cFuymKxSJJ/DoG2nutAZWRk+D/jANARCN0A0MX5foH3jV9s/ljoggULNGXKFEmNk0AlJyf7Q2phYaFGjBjhn0W5sLCwxQRdOTk5LSac2ttEXoWFhf7JtCZMmBDwS67NZpPBYGjzeNG28P0i3fRR+/29373VJzWO2W16PVr7Jbu1Cbia3oPk5GRNmTLFHwQyMjJkMBgkyT8e1ffapz3q2pvW7p/U8h4WFRX5a2j+HvZn0aJFMpvN/nBptVr3eZ/3dr32d63a+llsOqbXd02DndzNN+65rXzvp7S0NOCR9KysLGVkZCgrK0tms1lFRUX+kNrZHA6HbDZbiz8K+B6T39s1q6qqktT4SH1TvvOUlpa2aZ+2nktqnO+ipKQkYB9fmy98t/VcPr6fD/saHy7J//OlO0wQCKB7InQDQBfmC0oWi0VLly5Vbm7uPsc5FhUVac6cOTKZTCosLFRGRobMZrMKCgo0ZcoUzZw586Br8c3qnZWVpZKSEs2YMcPfmyc1/iKcnZ2ttLS0oN6ztCds+3q+mk9s1dr73V99+fn5ysrK8l/LGTNmtGmsuG+sqclk8k9wZ7FYtGjRIknSwoUL/WGhoKDAP9lWW6/bwdblM2PGjFZ7+fLy8mSxWGQ2m+VwODRlyhTNmDFDZWVl/vfgCzH7k5+fH9Bb6/v87W0yrL1dr/1dq7YqLCxUSkqKCgoKVFZWJrPZrMmTJx90sPVdO1+4a43vDwZNZy93OBwqKSnxhzaHw6H8/Hz/HwWmTJni/wp2HPzB8N1fk8nUoi0lJUVlZWWtHucLtc0/H77ru23btjbt09ZzSY2f18WLFysnJ0elpaX+n1+5ubn+UN3Wc/m22Ww2FRQUKCcnx/+1N2azudUx3wDQLkL9fDsA9HT7GtNtNpu92dnZAdvKysoCxpj6xmdardaA/UwmU4txir7xp03HSGZnZ7cY/9h8THF1dbVXkjcvL6/F+ZqPwzxQvnGbzb9MJpM3PT3dW11dHbB/a++3LfX5ztdUQUFBi+vRfDyo2WxucVxzvn+/+Tjl9qxrX8xmc4v7IMmbm5vr9Xr33M+D4bs/TT+fvvfV2nXZ3/Xa27Xyetv2WdzXOZte5wMZ1+v7HmyN7/3n5uZ6y8rKvGVlZa1+H7WHzMxMr9VqbfPX/saK+2pv7Vq39plp3t78Pvquk++4tuxzIPuVlJR4TSaT/2dA8599bT1XXl5ei5+HvnvW/OeJj8Vi2e/3OQAcrMgOzvQAgIPkG2PYfPyi2Wz29xw27X1s+oirzWaTw+Fo0UM8ZcqUg1r7uri4WFLj45rN69lX7+CByM3N9Y81lVrvnWuq6fvdX32+69G8bX9sNptsNttBL5HVUXU1l56ervz8fH+dvh5o3/X0PX0wYcIEZWVlyWq1tnk27by8PP+M3U17ki0WS4ue7mCv18HyfVb21nO7P23pIW/6eH1ubq7y8/OVk5PTrmPbQ3XdWuOb5TwnJ0czZsxQcXGxvz7f0ltt2aet+/l6tgsKCpSenh7wtEDTHui2nCszMzPgZ4m052dfcXFxqytDpKSkMKYbQIfh8XIA6KJ8ga35+EWpMQA0/wWxaYjqqGWKqqur5fV6A76aj8M8WL5xpr6v/Wntve2tvoO9Hu11Hdu7ruaysrICHjFftGiRrFar/zr6lk/y7TtixAhNmTKlTWFz8eLFcjgc/vHZvi/f+Nemwbszl8cqKipSRkaGfyx8Z8vNzVVhYWGXDmq++9Dafd5f3ZmZmf4/LmRkZKikpEQLFy6UtOcPbW3Zp637ZWRkKDs72/+HRLPZrKVLl6qoqChg/oC2/pvN+X6O7u19t+VnDgAcLEI3AHRRvl+YW/slsfnkSM1/YdzfL5jtWUsoNH+/+6vvYK9HsO+7o+pq7d+xWCz+Xldfr2FTFotFJSUlqq6uVl5enoqLizVv3rx9ntc3AVhZWVmLPxr4xmI37Z3trM9JRkZGQC/owYwLb+pgAldmZqbMZnPQTyl0JN8fsFasWBGw3Xd/9jdxXHZ2tqqrq/3rZvv+ENh03oa27LO//Xz1TJw4sUX9klr8Ya+t/2Zr73lvfxByOByt/oETANoDoRsAuihfT2XzR05LS0tVWlraYhmgpiwWi0wmU6sTkLWmeU9Y81l8faGutZDWFZbz2l99vuvR/FHg/U0k5nukuLXHfpu+b184aH4tOqqu1syYMUOLFy/29zzvbfkrk8mkzMxMWa3W/c7W7JsErbWgYjKZZLVaA2bmbsv12tu1ar6fT/MaHQ6HCgsLtXDhQn/wDZYvbB3oZzk3N9c/O3l7yMrKCph8bX9frU1k19z06dNbzDTvO661x6z3JTc3V9nZ2fv8I0Vb9mm+n+8eNp/IrK1/HGj+b7a2/6JFi2QymfYazKuqqujtBtBxOn8YOQCgqX1NpOabUCszM9O7ZMkSb15entdkMgVMErS3CaN8E45lZmZ6S0pK/P+Omk0A5du+ZMkSb3V1tbe
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1000x800 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Overall Torque R² Score: 0.999895\n",
|
||
|
|
"Overall Force R² Score: 0.999967\n",
|
||
|
|
"\n",
|
||
|
|
"Torque prediction is excellent!\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"y_pred_full = model_best.predict(X_poly_best)\n",
|
||
|
|
"torque_actual_full = y.iloc[:, 1].values\n",
|
||
|
|
"torque_pred_full = y_pred_full[:, 1]\n",
|
||
|
|
"\n",
|
||
|
|
"torque_r2_overall = r2_score(torque_actual_full, torque_pred_full)\n",
|
||
|
|
"\n",
|
||
|
|
"fig, ax = plt.subplots(1, 1, figsize=(10, 8))\n",
|
||
|
|
"\n",
|
||
|
|
"ax.scatter(torque_actual_full, torque_pred_full, alpha=0.4, s=15, color='#1f77b4')\n",
|
||
|
|
"ax.plot([torque_actual_full.min(), torque_actual_full.max()], \n",
|
||
|
|
" [torque_actual_full.min(), torque_actual_full.max()], \n",
|
||
|
|
" 'r--', linewidth=2, label='Perfect Fit')\n",
|
||
|
|
"\n",
|
||
|
|
"ax.set_xlabel('Actual Torque (mN·m)', fontsize=13)\n",
|
||
|
|
"ax.set_ylabel('Predicted Torque (mN·m)', fontsize=13)\n",
|
||
|
|
"ax.set_title(f'Torque: Predicted vs Actual (R² = {torque_r2_overall:.6f})', \n",
|
||
|
|
" fontsize=14, fontweight='bold')\n",
|
||
|
|
"ax.legend(fontsize=12)\n",
|
||
|
|
"ax.grid(True, alpha=0.3)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Overall Torque R² Score: {torque_r2_overall:.6f}\")\n",
|
||
|
|
"print(f\"Overall Force R² Score: {r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]):.6f}\")\n",
|
||
|
|
"print(f\"\\nTorque prediction is {'excellent' if torque_r2_overall > 0.99 else 'good' if torque_r2_overall > 0.95 else 'moderate'}!\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "1dd6a624",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Torque vs Roll Angle (Fixed Currents, Multiple Gap Heights)\n",
|
||
|
|
"\n",
|
||
|
|
"Examine how torque varies with roll angle for different gap heights. Roll angle is the primary driver of torque asymmetry."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "70b9b4e1",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAS1CAYAAADjrFn4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Xt4XFd1N/7v0d3WbSRfYgcnwaNA3JCAopGsOKEFailNItq0RSP3Ql+7SSUBDr1Q0MRQ2qa0GAlaWrBbJLt9lffl9xZLakkAK6aalAIl2JY0CBLAlGhMSHAcO5ZGN1vWZc7vj+GMNLqORjOz9pn9/TxPHvBctNY+a+nsOdpzzjFM0zRBRERERERERERERESkuDTpBIiIiIiIiIiIiIiIiKLBRQ0iIiIiIiIiIiIiIrIFLmoQEREREREREREREZEtcFGDiIiIiIiIiIiIiIhsgYsaRERERERERERERERkC1zUICIiIiIiIiIiIiIiW+CiBhERERERERERERER2QIXNYiIiIiIiIiIiIiIyBa4qEFERERERERERERERLbARQ0iIiIiIiICAHi9Xvh8Puk0iNaFfUxERJTauKhBRESUArxeL4qKitb0XyAQkE7bdgKBAAzDWPRfUVERqqur0dXVFfPPbmtrg2EY8Pv9Kz4Wy890uVwx5xVrzFhzlrBcXa1t19bWtux729ra4HK5lvzjWVdXF0pKSsI9YmlsbERRUREMw0BjY2NCxiTN2i7W2F0uFxobG8X6wu/3wzAMuN3uZV9j9e5af49X6oFo3msYhjL7Y6/Xi8bGRjidzvBj6/n9oDmrzR+x/gHejvvcpcTj92j+NnA6ndi7d6/ttwsREREtLUM6ASIiIlq/8vJyNDc3L3q8sbERDodjyeccDkcSMktNVVVVaG1tDf/b7/ejs7MTbrcbtbW16OzsFMxujpWHz+eD3++P+EMlLbawrj6fDz09PWhsbERnZyd6enoWvWdwcBA+nw9DQ0MRj3d1dcHtdqO5uRllZWXhP6xVV1djaGgIx44dS9nfwerqani9XtTW1oYXbfr7+9HR0YGOjg4MDw8LZxhfy/WA3fj9frjdbnR2di7Zm7H8fqSaQCAAj8cDr9cb3qe2traiqqoq6p8xfzsGAgH4/X60trbC5XKhubkZTU1NiUpfafH+PXI6nTh27BhcLhfOnz+fsvtbIiIiXXFRg4iIKAU4HA40NDQsetz6xu1Sz1HsnE5nxAKB0+lEVVUV3G43qqur0dLSIv6HqUAgAK/Xi+bmZng8HnR1dYnnlAxdXV1wOp0oKytb83uXqmttbW24ri6XC/39/RHvaW5uXnLR8PDhw2hoaIjY5n6/H16vF/39/THllyzr2YYlJSXhRb7a2tqI55qbm1PyW/3L9UC8racu0WhsbERVVdWyf6CP5fcjlfj9flRXV4e/KOB0OuH3+1FcXLymn7NwO5aVlYUXAD0eD6qqqpTePyRKIn6PamtrcfjwYXg8nogFOSIiIrI/Xn6KiIiIKE6qqqrgdDpx4sQJ6VTCfzxuamqCw+HQ5g869fX1cd/+1jerfT5f1Jcm8vl8i74ZbF1WRfVvDMe6DVtaWsLfOl+4oAGExq3DwlqiJKK3LdaCWyyXQ4vl98OOqqur4XQ60d/fj9ra2vBiRLwWIKx9tC776mQ5dOgQ2tralLnEGxEREcUHFzWIiIiI4kyFP56cOHEi/Ifluro6+P1+Xlt8HRoaGuBwOODxeKRTUdbhw4eXPWuM1GZ9Q34tl1GaL9V/P9ra2sJnICWSw+HgfjrOrHkwFc8SIyIi0hkXNYiIiDQWCATQ2NiIkpISFBUVwe12L/qDSldXV/gmx11dXeEb/lqsS3JYNzv1er1wuVwRN6de7kamC19n5eR2u1FUVISSkpKo/0jm8XiWveGudWmU+fnOv9FtS0tLVDFW09XVBb/fv+jbztFs53jy+/3w+XzYt28fAIRvjrzUN4Ct+i7McalvbPt8vvDNn60bBFu1qq6uXjGnWOtqvdfj8YRvul1SUrKoZm63O1z/lpaWcH3j+Ycsa3Fofo8t7G0r9vz/X1JSEvFzrHHMH8Nq22e138No379SndezDb1eLwKBAA4dOrTqa+eLprbWz3e5XPD7/eHXW79Libbatl1u/7aW35ehoaGY6xKPfZp1D5T1WOr3A4judz+abbXe34FoX7OU1tZWlJWVYWhoCG63GyUlJSgpKYnr/iUQCCAQCCy6Sfta546V5sL5c+5a9v1r+axgvbaoqCji582vcUlJyaKzepb6PYp2/7CasrIyLe75QkREpBMuahAREWkqEAhg586d4fsuHDt2DIFAAC6XK3yZHCD0xzbrj2lutxvFxcXhPzL5/X6UlJTA6/Xi0KFDcLvd8Hg8Ee9fC7/fj507d8Lv9+PYsWPweDxoa2uL6g9P1h/wOzo6Fo3T6/WGn7f+oNPT04POzk5UVVWt+Y8d1qKB9Z91U2i3242qqqqIS+xEu53jyfpjkfWta+t/l7o0jFVfl8sVvlZ8VVUV2traFv3R3eVyoby8HP39/aiqqkJjYyMqKirwzDPPrHgt9PXUFQjV1OfzwePxhG9M7PF4It7f3Nwcvp5/bW0t+vv70d/fj7q6uqhiRMNanOjr61v2NQ0NDYvyWNhfra2t6O/vD5/REM32We33MNr3r1Tn9WxDq5fXeimeaGprjdH6oygQ+sNteXl5+HcvWoFAIOJ3d/5/S90PItbeXevvy3rqst59mnWz6oqKiqjfs5Slfj+i2X7Rbqv1/g6sZz9kLdZUV1ejuro6fE+NxsZGeL3edW03i5WHtQgQ69xhvf/w4cOLxuDz+cLPr2Xfv5bPCtbPO3bsWPjnVVdXw+1249ChQ+HFdbfbvepZjdHuH1ZTXl6+4n6biIiIbMgkIiKilAXALCsrW/K52tpa0+l0Lnq8rKws4j2tra0mABOA2dPTE/HaqqoqE4A5ODgYfmxwcHBRXOtnzH/dUrGqqqpMh8MR8Zqenh4z2o8sTqfTrKqqinhsfuz+/v4l84jW8PBweFss/M/pdJqdnZ2L3rPW7Tw/t+W222qcTueiulu16u/vj3jcitHQ0BDxOICIbdnU1LSoNk6nc9H7lsp5vXVdSm1t7ZLvdzgcZlNT05p+llXXhWNZqLOz0wRgtra2hh9brkZL/Tzr/QtrEM32We33MNr3r1Zn04xtGzY0NKzrd2u+pWpr5b/wd2ypfdBSrP1SNP/Nj7GWbTs/h7X+vsRal/Xu0+b/jKX2X6a5vt+PaLbfWrdVrL8Dse6H5u/3F27npfa1q/2cqqoqs7+/P/xfZ2dnuI+bm5vDr1/P3LHUWJuamkwA5vDwcMT7Vuu9teax8Oc5HI5F+z1ru0ezL11opf3Dcu9tbm6OGDsRERHZH8/UICIi0lAgEEBXV9eS33Zsbm4Of3N5vqampkXXW/d6vWhoaIi4XIbT6Yz491pysn6edRmOQCCA8vJyOByOqL4NW1tbG74MjsW6bIjT6QzfoLm5uXld971oaGiAaZrh/6yfvfDyLbFs5/Wy7p1hnZliWekSVPOft1iXWrH4fD6Ul5dHvMbpdK56KZR41HUp1rfK11PHtbK2Ryz9vZy1bp+Fv4drff9qdY6V9bu1VD28Xm/Ef6tZqbYLzwSxvsUf7Q2qa2trI3535/+38HdjPb271t+XWOsSj31avPp64c+JdvutdVvF8jsQj/1QVVXVom1UW1sLn8+3pm0//xKNLpcL9fX1AEJn2lhn+a137mhsbAyP2dLV1YXa2tpwz1hW6r1Y8lj486z5cf7v7vweWatY9v3WmOOxryMiIiI1cFGDiIhIQ9ZlGBb+IWn+Ywsv1bDwGvDWHzIW3hMDwKI/mkRj/j0JrGtxW/9Zl0dZjXVZDesSVNalZqzHnU4nGhoa0NbWFtO9HZZj/XFn4fXVY9nO62X9Ydb6o6D1n/VHpIWX57Ks9gfNpf7IOP/nLicedQXmrq3ucrlQVFQkckNi6xJA8VzUWOv2Wfh7uNb3xzP3+aw/NC7s5/mX7LH+W5jTempr/aF0cHBwnSNYbD29u9bfl1jrEo9
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1600x1200 with 4 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Showing torque variation with roll angle at 4 different gap heights.\n",
|
||
|
|
"Gap height affects the magnitude of torque but the relationship with roll angle remains consistent.\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was written by AI - see [9]\n",
|
||
|
|
"# Select 4 different gap heights to visualize\n",
|
||
|
|
"# Using currL = -15A, currR = -15A configuration\n",
|
||
|
|
"gap_heights_to_plot = [8, 10, 15, 21] # mm\n",
|
||
|
|
"\n",
|
||
|
|
"fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n",
|
||
|
|
"axes = axes.flatten()\n",
|
||
|
|
"\n",
|
||
|
|
"for idx, gap_height in enumerate(gap_heights_to_plot):\n",
|
||
|
|
" ax = axes[idx]\n",
|
||
|
|
" \n",
|
||
|
|
" # Fixed current configuration\n",
|
||
|
|
" currL = -15\n",
|
||
|
|
" currR = -15\n",
|
||
|
|
" \n",
|
||
|
|
" # Get actual data for this gap height\n",
|
||
|
|
" gap_data = magDf[(magDf['GapHeight [mm]'] == gap_height) & \n",
|
||
|
|
" (magDf['currL [A]'] == currL) & \n",
|
||
|
|
" (magDf['currR [A]'] == currR)]\n",
|
||
|
|
" \n",
|
||
|
|
" # Create fine grid for smooth curve across roll angles\n",
|
||
|
|
" roll_fine = np.linspace(-5, 5.0, 500)\n",
|
||
|
|
" X_condition_fine = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [currL] * 500,\n",
|
||
|
|
" 'currR [A]': [currR] * 500,\n",
|
||
|
|
" 'rollDeg [deg]': roll_fine,\n",
|
||
|
|
" 'invGap': [1/gap_height] * 500\n",
|
||
|
|
" })\n",
|
||
|
|
" \n",
|
||
|
|
" # Predict with best degree polynomial\n",
|
||
|
|
" X_condition_poly = poly_best.transform(X_condition_fine)\n",
|
||
|
|
" y_condition_pred = model_best.predict(X_condition_poly)\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot polynomial curve for TORQUE (index 1)\n",
|
||
|
|
" ax.plot(roll_fine, y_condition_pred[:, 1], '-', linewidth=2.5, \n",
|
||
|
|
" color='#ff7f0e', alpha=0.8, label=f'Degree {best_degree} Polynomial')\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot actual torque data\n",
|
||
|
|
" ax.scatter(gap_data['rollDeg [deg]'], \n",
|
||
|
|
" gap_data['YokeTorque.Torque [mNewtonMeter]'],\n",
|
||
|
|
" s=100, marker='o', color='#9467bd', \n",
|
||
|
|
" edgecolors='black', linewidths=1.5,\n",
|
||
|
|
" label='Ansys Data', zorder=5)\n",
|
||
|
|
" \n",
|
||
|
|
" # Formatting\n",
|
||
|
|
" ax.set_xlabel('Roll Angle (deg)', fontsize=11)\n",
|
||
|
|
" ax.set_ylabel('Torque (mN·m)', fontsize=11)\n",
|
||
|
|
" ax.set_title(f'Gap Height = {gap_height} mm (currL={currL:.0f}A, currR={currR:.0f}A)', \n",
|
||
|
|
" fontsize=12, fontweight='bold')\n",
|
||
|
|
" ax.legend(fontsize=10, loc='best')\n",
|
||
|
|
" ax.grid(True, alpha=0.3)\n",
|
||
|
|
" \n",
|
||
|
|
" # Add vertical lines at actual roll angle data points\n",
|
||
|
|
" for roll_point in gap_data['rollDeg [deg]'].unique():\n",
|
||
|
|
" ax.axvline(roll_point, color='gray', linestyle=':', alpha=0.2, linewidth=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.suptitle(f'Torque vs Roll Angle at Different Gap Heights (Degree {best_degree} Polynomial)', \n",
|
||
|
|
" fontsize=15, fontweight='bold', y=1.00)\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Showing torque variation with roll angle at {len(gap_heights_to_plot)} different gap heights.\")\n",
|
||
|
|
"print(\"Gap height affects the magnitude of torque but the relationship with roll angle remains consistent.\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "02d02793",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Torque vs Current Imbalance (Fixed Roll Angle, Multiple Gap Heights)\n",
|
||
|
|
"\n",
|
||
|
|
"Now examine how torque varies with current imbalance at different gap heights for a fixed roll angle."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "3b8d7f6a",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAS1CAYAAADjrFn4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs/XtYW9edL/6/t7gZMFhg40tiJ0Hk4iRO6wiwc+k1Fp2mbqdnWuHMzJmOM0mB03HOXPKdgeNzzvP9Pvl+n994YHpm5swkZw64nTpzejPQJmkbNy3KtE3bNDagkMZJnAty7k7sGGTAXMXevz92toLQBQHSZ21J79fz8CTW3pIWby0ttLT2WkszDMMAERERERERERERERGRzTlUF4CIiIiIiIiIiIiIiCgZHNQgIiIiIiIiIiIiIqKMwEENIiIiIiIiIiIiIiLKCBzUICIiIiIiIiIiIiKijMBBDSIiIiIiIiIiIiIiyggc1CAiIiIiIiIiIiIioozAQQ0iIiIiIiIiIiIiIsoIHNQgIiIiIiIiIiIiIqKMwEENIiIiIiIiIiIiIiLKCBzUICIiIiIiomXz+Xzw+/2qi0GUs/geJCKiXMVBDSIiyno+nw8VFRXL+gkGg6qLnfG6urpQW1sLTdOgaRpqamrQ2NgIn8+numjiAoEANE1DY2NjSh+3q6sLmqYhEAjY8vGyUTAYDNfpxT+1tbXo6uqKe1/rPRHrC6je3l7U1NRA0zRUVFSEb29paUFFRQU0TUNLS0tafifVFrYVFRUVqK2tRUtLi7J6mMz71Xqv9Pb2LuuxE9WBZO6raZpt/j75fD60tLTA5XKpLkqEdLW32cCqQ7W1tSl93HjtYkVFBRoaGpb9Plko1t8l/q36gMvlwp49e5gFERHlnHzVBSAiIkq3uro6tLe3R93e0tICp9MZ85jT6RQoWXYKBoPhwYvm5mYcPHgQIyMjGB4eRldXFyorK+HxeFQXk2hVPB4POjs7w//2+/3o6+tDS0sLenp60NfXF3Wf4eFh+P1+jIyMRNze29uLxsZGtLe3w+12h7+camhowMjICA4fPpy1bVJDQwN8Ph+8Xm940GZwcBDd3d3o7u7G6Oio4hKmVrw6kGkCgQAaGxvR09OTtXVTpWAwCJ/Ph7a2NgwODkZlbLUZsXR2dqK5uTnmsZ6eHgBmexUIBFI+ILW4XQwEAujp6UFjYyO8Xm/4+bNNS0sLuru7EQwG4fV6k26zk30dA4FAeKB3ccYulwuHDx9GbW0tTp8+zfcjERHlDA5qEBFR1nM6nTE7+NYVpvE6/7QytbW1CAQC6Ovrixq8aG9vz+irCXt7e+FyueB2u1UXhVJgNa+ny+WK+ELQ5XLB6/WisbERDQ0NqK2txeDgYMR92tvbYw6iHjp0CM3NzWhtbQ3fFggE4PP5MDg4aOv6tpoMa2pqwl96er3eiGPt7e0JZ71kqnh1INXS3Va1tLTA4/FwgDrFAoEAampqAJifXZaalRNrUCneQIU1UNLe3o62tjb09vZGtDmpEKtd9Hg84Xaxo6Mj5c+pWm1tLYLBIA4ePAiXy4VDhw6hurp6WQOyS72ObW1t6OzshMvlCr92C9tMr9eLQ4cOhc8jIiLKBVx+ioiIiFKmra0t7oCGxW5LlSxHU1MTjh49qroYlCLpeD2tq2j9fn/SS674/f6oL7Ss5YnsftXtSjPs6OhAIBBAZ2dn1IAGYP7e2fblp6R0tlXWgFu2LoemksvlwvDwMAzDSOqCC2tgaeFPvL+x1iBha2srnE6n6JffVrmy7e9nV1dXeJZea2srvF4vHn/8cQSDQbS1tSX9OEu9jgsHt+INdB08eBBdXV22WZ6OiIgo3TioQURERCkRDAbR0dERvjKTKFc1NzfD6XQu60utXHPo0KG4s+jI3qyZJmzn0yNdA/9Hjx4NDyDu27cPgUBAfOZktn3h3tPTEzUAYbVrq9lHZLH29nY0NjaGZ/HEGgi2bsvGGW5ERESxcFCDiIgojmAwiJaWFtTU1KCiogKNjY1RXwD09vaGN/Xt7e0Nb3BrCQQCaGhoCG+W6fP5UFtbG7FJZ7wNLxefZ5WpsbERFRUVqKmpSfpL07a2trgbzFpL5Sws78KNjzs6OpJ6ju7u7vBzJSvZ3z1Rzku9BktlZt1/8eu98DEaGxvD+XV0dITzWe2XB4uf29qo3npuv98f3kC5pqYm7pck1nrbscpuZdDW1hbeiLqmpibp1zXZ+yaT40LWOfEec6V1PZnypuv1XMj60nDhe25xfbeee+H/W19aWazfY+HvkGydtv7fbu8Jn88XXq5lOZKti1Y7GwgEwudbbXi6LZVtvDZv4Xvd2nDeepyGhoaIc0dGRlb8uqymjbdYe6Akkuj9LdHu55JgMLjk4EQgEIDf78cdd9wBAOH3gtRsjd7e3vDfqYWS+ZxlZz6fL+YSb1b7s5xBnESvo9vtxuDgIIaHhxO+Zm63O+Z+TkRERNmIgxpEREQxBINBVFdXh9efPnz4MILBIGpra8PLwgDml0vWl0eNjY2orKwMfwFlrY3t8/lw8OBBNDY2oq2tLeL+yxEIBFBdXY1AIIDDhw+jra0NXV1dSX3Za32RYQ08LPw9fT5f+Lj1hVJfX1/4CsRkO8jDw8MA0nOVaaKcl3oNlsrMun9tbW1443iPx4Ourq7wF3Ht7e3h/RG8Xi8GBwcxODiIffv2peT3sp778OHD4eduaGhAY2MjDh48GP4So7GxMeaXJI2NjVFlX/gFbnd3N/x+P9ra2sKbWbe1tSVVd5K9bzI5Amadq6mpQXd3d3hTbWs9cOu9sZq6nkx50/V6LmQNTgwMDMQ9p7m5Oaoci99vnZ2dGBwcjNgwNtk6bdf3hPU6L3e/h2TrovUFrtWetbW1oa6uLuGmvLEEg0H4/f6YP4v3S7GedyX11sq6rq4Og4OD8Hg8aGlpQX19PR5//PGoPThW87qspo23yhoIBFBfXx/3+FLv72SttN3PJQsH0CoqKuDz+WKeZw2IW7NrrP+mcjYB8MF7z/qx3nONjY3weDwRS8ol+znLrqy/xevXr486VllZCQBJD9Ak+zoupa6uLuHfHCIioqxiEBER5SgAhtvtjnnM6/UaLpcr6na32x1xn87OTgOAAcDo6+uLONfj8RgAjOHh4fBtw8PDUc9rPcbC82I9l8fjMZxOZ8Q5fX19RrJ/zl0ul+HxeCJuW/jcg4ODMcuRrObm5mXfP9nfPVHOS70GS2Vm3b+5uTniPABReTmdTqO1tTXp389ive5erzdm2Rc/t9PpNAAYg4ODUeXu7OyMuv/i3zuZ18Lr9UbVnXivx3Luu1SOzc3NhtPpNEZHR+M+/mrrejLlNYyVvZ6jo6Mxf8/Fenp64r5ei/ON9XjW/RfWAcNYXp2263tiJW1FPInqYk9PT8TtsdrkWKz3azI/C59jOdkuLENra2vU/VwuV1T+q31dVtvGL3yMxdlaknl/p7vdj9feLldzc3O4TMn+LHy/r1Zra6sBIGaWVvvgcrmMwcFBY3R0NPxeWNxmGIZZnxZ/3rHeD7HOXy6rXYz143K5YtaX5X7OWlhfkv1bZUnHa2nVs1jnWe/7xXVzseW+jktpb2+PW2eIiIiyDWdqEBERLRIMBtHb2xvz6tr29vbwFYgLtba2Rq0v7vP50NzcHDFzweVyrWgmgzWjorm5GcFgMPxTV1cHp9OZ1FV9Xq83vOyLpbOzE263Gy6XK7whcXt7+4rWvbbun86lI2LlHO/YcjNbfAW32+3GyMhI6n+JGBY/t/V6LLyS3ao3sV6bxXXKqruJrsK1rrReyWud6L5L5djV1RXecyKWVNT15ZQ3XazfOZUzl5abjV3fE9ZrH+v18Pl8ET9LSfTaLp4JYs14SPbqdK/XC8MwYv4sXgJmNfXW7/ejrq4u4jaXyxW3LV3p67LaNh5Yul4v9f5eieW0+6lkzZJazo/UHjEulwterxd9fX1wu91wOp3o6ekBEL0EpLV3hjUj05KOJaiam5sj3ifW37LFy5Wt5HPWatj1tVzO65gM630n9dmFiIhIJQ5qEBERLWJN3V/8JdP
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1600x1200 with 4 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Showing torque variation with current imbalance at 4 different gap heights.\n",
|
||
|
|
"All data uses currL = -15A with varying currR to create different imbalances.\n",
|
||
|
|
"Current imbalance (currR - currL) affects torque magnitude and direction.\n",
|
||
|
|
"Positive imbalance (currR > currL) and negative imbalance (currR < currL) produce opposite torques.\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was written by AI - see [10]\n",
|
||
|
|
"# Select 4 different gap heights to visualize\n",
|
||
|
|
"# Using roll = 0.5 degrees (small but non-zero torque)\n",
|
||
|
|
"gap_heights_to_plot = [8, 10, 15, 21] # mm\n",
|
||
|
|
"roll_target = 0.5\n",
|
||
|
|
"\n",
|
||
|
|
"# Use a fixed reference current for currL\n",
|
||
|
|
"currL_ref = -15\n",
|
||
|
|
"\n",
|
||
|
|
"fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n",
|
||
|
|
"axes = axes.flatten()\n",
|
||
|
|
"\n",
|
||
|
|
"for idx, gap_height in enumerate(gap_heights_to_plot):\n",
|
||
|
|
" ax = axes[idx]\n",
|
||
|
|
" \n",
|
||
|
|
" # Get actual data for this gap height, roll angle, AND fixed currL\n",
|
||
|
|
" gap_roll_data = magDf[(magDf['GapHeight [mm]'] == gap_height) & \n",
|
||
|
|
" (magDf['rollDeg [deg]'] == roll_target) &\n",
|
||
|
|
" (magDf['currL [A]'] == currL_ref)]\n",
|
||
|
|
" \n",
|
||
|
|
" # Calculate imbalance for actual data points (all have same currL now)\n",
|
||
|
|
" actual_imbalances = []\n",
|
||
|
|
" actual_torques = []\n",
|
||
|
|
" for _, row in gap_roll_data.iterrows():\n",
|
||
|
|
" imbalance = row['currR [A]'] - row['currL [A]']\n",
|
||
|
|
" torque = row['YokeTorque.Torque [mNewtonMeter]']\n",
|
||
|
|
" actual_imbalances.append(imbalance)\n",
|
||
|
|
" actual_torques.append(torque)\n",
|
||
|
|
" \n",
|
||
|
|
" # Create fine grid for smooth curve across current imbalances\n",
|
||
|
|
" # Range from -10A to +10A imbalance (currR - currL)\n",
|
||
|
|
" imbalance_fine = np.linspace(0, 35, 500)\n",
|
||
|
|
" \n",
|
||
|
|
" # Vary currR while keeping currL fixed\n",
|
||
|
|
" currR_fine = currL_ref + imbalance_fine\n",
|
||
|
|
" \n",
|
||
|
|
" X_condition_fine = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [currL_ref] * 500,\n",
|
||
|
|
" 'currR [A]': currR_fine,\n",
|
||
|
|
" 'rollDeg [deg]': [roll_target] * 500,\n",
|
||
|
|
" 'invGap': [1/gap_height] * 500\n",
|
||
|
|
" })\n",
|
||
|
|
" \n",
|
||
|
|
" # Predict with best degree polynomial\n",
|
||
|
|
" X_condition_poly = poly_best.transform(X_condition_fine)\n",
|
||
|
|
" y_condition_pred = model_best.predict(X_condition_poly)\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot polynomial curve for TORQUE (index 1)\n",
|
||
|
|
" ax.plot(imbalance_fine, y_condition_pred[:, 1], '-', linewidth=2.5, \n",
|
||
|
|
" color='#ff7f0e', alpha=0.8, label=f'Degree {best_degree} Polynomial')\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot actual torque data (now only one point per imbalance value)\n",
|
||
|
|
" ax.scatter(actual_imbalances, actual_torques,\n",
|
||
|
|
" s=100, marker='o', color='#9467bd', \n",
|
||
|
|
" edgecolors='black', linewidths=1.5,\n",
|
||
|
|
" label='Ansys Data', zorder=5)\n",
|
||
|
|
" \n",
|
||
|
|
" # Formatting\n",
|
||
|
|
" ax.set_xlabel('Current Imbalance (currR - currL) [A]', fontsize=11)\n",
|
||
|
|
" ax.set_ylabel('Torque (mN·m)', fontsize=11)\n",
|
||
|
|
" ax.set_title(f'Gap Height = {gap_height} mm (currL = {currL_ref}A, roll = {roll_target}°)', \n",
|
||
|
|
" fontsize=12, fontweight='bold')\n",
|
||
|
|
" ax.legend(fontsize=10, loc='best')\n",
|
||
|
|
" ax.grid(True, alpha=0.3)\n",
|
||
|
|
" ax.axhline(0, color='black', linestyle='-', linewidth=0.5, alpha=0.3)\n",
|
||
|
|
" ax.axvline(0, color='black', linestyle='-', linewidth=0.5, alpha=0.3)\n",
|
||
|
|
" \n",
|
||
|
|
" # Add vertical lines at actual imbalance data points\n",
|
||
|
|
" for imb in set(actual_imbalances):\n",
|
||
|
|
" ax.axvline(imb, color='gray', linestyle=':', alpha=0.2, linewidth=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.suptitle(f'Torque vs Current Imbalance at Different Gap Heights (currL = {currL_ref}A, Roll = {roll_target}°)', \n",
|
||
|
|
" fontsize=15, fontweight='bold', y=1.00)\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Showing torque variation with current imbalance at {len(gap_heights_to_plot)} different gap heights.\")\n",
|
||
|
|
"print(f\"All data uses currL = {currL_ref}A with varying currR to create different imbalances.\")\n",
|
||
|
|
"print(f\"Current imbalance (currR - currL) affects torque magnitude and direction.\")\n",
|
||
|
|
"print(f\"Positive imbalance (currR > currL) and negative imbalance (currR < currL) produce opposite torques.\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "d71ac3cc",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Torque Segment-Wise Performance\n",
|
||
|
|
"\n",
|
||
|
|
"Calculate R² scores for torque predictions across all segments, similar to what we did for force."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "902343b2",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Torque Prediction Analysis (across 637 segments):\n",
|
||
|
|
"\n",
|
||
|
|
"Best R² Score: 0.999964\n",
|
||
|
|
" Parameters: currL=15A, currR=-15A, roll=-1.5°\n",
|
||
|
|
"\n",
|
||
|
|
"Worst R² Score: -0.638045\n",
|
||
|
|
" Parameters: currL=-10A, currR=-10A, roll=0.0°\n",
|
||
|
|
"\n",
|
||
|
|
"Mean R²: 0.981564\n",
|
||
|
|
"Median R²: 0.999392\n",
|
||
|
|
"Std Dev: 0.124378\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAJNCAYAAAAs3xZxAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXShJREFUeJzt3X10G/d95/sPqAdKtkUOIEuxEys2B3brdJNuDFBJ03TbxgLdrHO3jdcA1bRN67s1gTinTbvqDRG26d361A0N1pu27p7YIJ093vasUxFoHrrnpqkBOclu0yYhASuN2zxiaMeO7VAWAJGyJVoUcf9QMCVEAARBDEiC79c5PBJmfpj5zgAD8oP5zW9cpVKpJAAAAAAA0HJdG10AAAAAAACditANAAAAAIBDCN0AAAAAADiE0A0AAAAAgEMI3QAAAAAAOITQDQAAAACAQwjdAAAAAAA4hNANAAAAAIBDCN0AANRhWZYGBgbkdrvldrsVCoVkWdZGlwUAALYIQjeAmrLZrFwu14oft9utSCSiYrFY9Xnj4+MKhUKO1XX58sfGxuR2ux1bX711b0bRaFRut1sul0tjY2M126XTaTtINvLTKUGz1vva7/crGo2ueF97vV75fD4lEgnFYjFls1n5/f6a7//LFYtFRSIReb1eez2RSKRj9udyfr9fLpdL2Wx2o0txVCPbefnn0lo+p+q9R+sd02jMdjomAWwOOze6AACbXzwe1+DgoCQpn88rnU4rFoupr69PJ06ckM/nq2ifyWSUTCYdq8fp5W/WdTciEoloenpaiURCHo+nbtv+/n4lEomKaYlEQuPj44rFYiteV9M0W17vRrr8fZ3NZjU6Oqrx8XFlMhl7exOJhILBoP28QCAgr9erdDpdMb2aYrGovr4+eTweRaNRmaYpy7IUi8WUz+dX7P+tzLIsZbNZGYaheDyueDy+0SU5op3bWe2zNxqN6vjx48pkMo6tt5Ntp2MSwCZSAoAaMplMSVIpkUhUnR8IBEqSSoVCYV3rSSQS61pGLBYrGYaxrhout96aNkq916sR8Xi8JKmUyWRaWNXmstr72ufzlXw+X83np1KphvdRMBgsmaZZdV4ul2us4C0iFouVAoFAKRwOt/x43Ewa3c7LP5fW8jlV7z2ay+VKkkrDw8NrL34Tavdn7XY6Jpfbqr/TgE5B93IATSufEYhGo+taTigU0vT0dCtKapnNWNNqyt2dDcPY0Dq2ukAgULPbcPlMYyKRWNEToFb7QCBQdV4n9hwIhUIKhUIqFotKp9MbXZIjNno7TdNUMBjsmP3b7s/a7XRMLrcVf6cBnYTQDaBphmEoHA5rfHy84etbgc0unU6v+OM7m81qYGBAlmUpk8ms2q28zDTNjglH9WSzWVmWpUAgYAeaTuymu5m2k8/c5myXYxLA5kLoBrAuAwMDklQxAM3lAwal02l74CG3262BgQFls1mFQiG5XC57OeXBgpYvp7z8SCQit9utYrFYc0CidDptL6d8ze1y5QHGLn+Oy+Wy62+kpsuXUQ5kbrdbXq93xUBHY2Nj8vv9drvywD1rHWyq3nqW11Vex/j4+JqWv94aynVUe80kKZlMyuv12iOAJ5PJin0vNfYalSWTSft95ff71/2HdPk9mc1mK67TjUQiCoVCikQi9h/s6XS6odBTHpyp0fqW79/lx0q1+c3s/1r7q9Yx2qjjx4/LNE37y4pAIFD3/VdvO+ttw2rbX2871ruNzWynE8pn18uh37IshUIhe/DEap8ttfZpo88tT7/8/VMsFu3ne73eqvui1vtutc/aes+tt02rWesxudrnzPj4uD0gWyQSUTQatds7uf8uX3at3y2r7edWHBcAGrDR/dsBbF6rXfu6vE08HrenLb92sVAolCSVYrFYKZfLlVKpVGl4eLiUSqVKhUKhYh2FQqHimrNYLFYyTbPk8/lKpmmWYrHYiuWXH0sqBQKBUiqVKqVSqZLP5ytJqrhGb3h4eMU1leXrc8vtGqlp+TLKzx8eHi5lMplSIpEomaZZCgaDFeuVVPL5fHZ9pmmu6brXRtaz/LVo9tq9etd0N1JDrdcskUhUvEbxeLxkGEZTr9HyOuPxeCmTydjvgUauySzvp2o/gUCgYtvL799qP8vf8/WUa1u+jlQqVbOuYDBov0+Gh4fta3fXs//r7a96x2ijDMOouMa4vL5qnx2rbWetbVht+1f7rFnvNq51O1t9TXehUKj47Cgf47FYrBQOh0upVKqUy+VKwWCwYn65Ta3P0tWeW+/zy+fzlRKJRCmVStljfDR6nK72WbvaMV5rmxrR6DG5Wg3lz7VUKlXKZDKlQCBQMk2zlMvl7DZO7b/Vlr38fVNrP7fquACwOkI3gJoaCd3lQX1qhe7yH8q1lH/pV/slvzxMXz69Wui+fLmGYZTC4bA9rdFAt1pNy5dhmuaKAY3K+6S838p/GC1fXvmPqUbDcSPrKdft1EBqjdRQ6zUr/4G5XHm/NPPFSLXQOzw8XPF611J+X5f/0MzlclVraaVCoVCKx+OlYDBYM7RfHqAv1+z+X21/rXaMrqa8P6t9WVFte1bbzlrvodW2v952rHcbS6W1b2crQvflP4ZhlILBYN3PjWqvd6192shz631+LX89Lv990chxWuuztpHnNrpN9ba13jHZSA3lAfUu3wfVvrRwYv81+rul1n5uxXEBoDF0LwewLuVuaLUGoOnv75d06b624+PjTd0HtZlb8hiGoUAg4Oi1e+XrOyORSMV00zTl8/l0/PjxiunLB+8p7698Pt/y9ThhrTUsf80sy1KxWFQsFqtoU+4aulblwYAikUjFPYzHxsbWNFBQuZuwaZqKxWIyDGPdgwLWUh7/IJFIqFAoyOfzVdzr3rKsqvu3bD37f7X9td5jNB6PyzAMmaapYrFob5PP51txe73VtrPWNjSy/fW2o1WfQ41uZ6vEYjEVCoWKn0QiUXewxPK8XC63Yt5qn6X1nlveh9K/fn4tP4bLAwuWP9PWc5yu5bnN3rJttWOykRoa+fwuc3L/Nfu7pRXHBYDGELoBrMvU1JSkyj8oljMMw76fbCQSkdfr1cDAwJoGAWp2RFnTNNf0R9Falf/wqXY/7PK9X7fSelpZw/LXrDyv1SMDFwoFlS712LJ/1nPv4lgspmQy6fj+NAxDExMTkv51v662j9az/8tq7a/1HqOTk5MqFov29dnln/IXcssD6VreC8vbNLL99bajFZ9Da9nOVjEMo+KnlnQ6rVAoZI+ZUEu1/d7oc6utv9rrcbn1HKeNPLcVnyvVjslGaohEIpqcnLSvzx4dHZXP56u6rzZi/62mFccFgMYQugE0rVgsanx8XMFgsO4fhD6fT5lMRoVCQfF4XNPT0xodHW1oHeu5/ZVlWY7eAqa87Goh7fJ1r2c71rIep6xnW8t/WLYqzNarZT3C4bBM02zoLOx6lfdR+Y/b1bZpPfu/kf3V7DFaDhu5XG5FMCgUCpIqz0Q2+tqtZRuWb3+97VjP59Bat7OdyrcvGxgYUCqVsuu5XLXPoPU8t950aX3HabPvk/VY6zEpXeoRYJqmPRhasVjUiRMnGq6zFftvvftgPccFgMYRugE0bWhoqGq34VrK3fmW3wf58j90WuXyEX6XT1+u2iitjdYUCARkGMaKP7az2ayy2ayOHj265ro3cj1O1VA+83P5+6TWH5OrvUblLsXV/jBc7/soFovZo5O3wvLuqsuV92P5/Vnu5l4tuBWLxXXt/7Xsr2rHaD3l+5VX++Jn+SUey4NMve2sZa3bX2871rqNzWxnuxSLRSWTSU1MTNhfGrXjuY1o5H1X67PWyWN8LcfkajVks1n7EoBSqaRUKtWyLwJauQ8a+Z3WzHEBoHE7N7oAAJtfPp+3f1nn83lls1mNjo7KsiylUqm6f6wlk0mNjo5qZGTE7gaaTqc1MjJitzEMw74VTzwet6+vXauBgQFFo1EVi0X72tzl6/F6vZIunbXq7+9XOp2u+Y1+ozVNTEwoFApJunTWyLIsRaNRBQKBhu/l3Ih2rcepGkZGRhSNRhWJRBSJRDQ9PV31jHKjr9H
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1000x600 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Following generated by AI - see [11]\n",
|
||
|
|
"# Calculate R² score for torque in each 13-row segment\n",
|
||
|
|
"torque_segment_scores = []\n",
|
||
|
|
"torque_segment_info = []\n",
|
||
|
|
"\n",
|
||
|
|
"for i in range(num_segments):\n",
|
||
|
|
" start_idx = i * 13\n",
|
||
|
|
" end_idx = start_idx + 13\n",
|
||
|
|
" \n",
|
||
|
|
" # Get actual and predicted torque values for this segment\n",
|
||
|
|
" segment_actual_torque = y.iloc[start_idx:end_idx, 1].values # Torque column\n",
|
||
|
|
" segment_pred_torque = model_best.predict(poly_best.transform(X.iloc[start_idx:end_idx]))[:, 1]\n",
|
||
|
|
" \n",
|
||
|
|
" # Calculate R² for this segment\n",
|
||
|
|
" segment_r2_torque = r2_score(segment_actual_torque, segment_pred_torque)\n",
|
||
|
|
" torque_segment_scores.append(segment_r2_torque)\n",
|
||
|
|
" \n",
|
||
|
|
" # Store segment info\n",
|
||
|
|
" segment_data = magDf.iloc[start_idx:end_idx]\n",
|
||
|
|
" currL = segment_data['currL [A]'].iloc[0]\n",
|
||
|
|
" currR = segment_data['currR [A]'].iloc[0]\n",
|
||
|
|
" roll = segment_data['rollDeg [deg]'].iloc[0]\n",
|
||
|
|
" torque_segment_info.append({\n",
|
||
|
|
" 'start': start_idx,\n",
|
||
|
|
" 'end': end_idx,\n",
|
||
|
|
" 'currL': currL,\n",
|
||
|
|
" 'currR': currR,\n",
|
||
|
|
" 'roll': roll,\n",
|
||
|
|
" 'r2': segment_r2_torque\n",
|
||
|
|
" })\n",
|
||
|
|
"\n",
|
||
|
|
"# Find worst and best performing segments for torque\n",
|
||
|
|
"worst_torque_idx = np.argmin(torque_segment_scores)\n",
|
||
|
|
"best_torque_idx = np.argmax(torque_segment_scores)\n",
|
||
|
|
"worst_torque_segment = torque_segment_info[worst_torque_idx]\n",
|
||
|
|
"best_torque_segment = torque_segment_info[best_torque_idx]\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Torque Prediction Analysis (across {num_segments} segments):\")\n",
|
||
|
|
"print(f\"\\nBest R² Score: {max(torque_segment_scores):.6f}\")\n",
|
||
|
|
"print(f\" Parameters: currL={best_torque_segment['currL']:.0f}A, currR={best_torque_segment['currR']:.0f}A, roll={best_torque_segment['roll']:.1f}°\")\n",
|
||
|
|
"print(f\"\\nWorst R² Score: {min(torque_segment_scores):.6f}\")\n",
|
||
|
|
"print(f\" Parameters: currL={worst_torque_segment['currL']:.0f}A, currR={worst_torque_segment['currR']:.0f}A, roll={worst_torque_segment['roll']:.1f}°\")\n",
|
||
|
|
"print(f\"\\nMean R²: {np.mean(torque_segment_scores):.6f}\")\n",
|
||
|
|
"print(f\"Median R²: {np.median(torque_segment_scores):.6f}\")\n",
|
||
|
|
"print(f\"Std Dev: {np.std(torque_segment_scores):.6f}\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Create histogram of torque R² scores\n",
|
||
|
|
"fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n",
|
||
|
|
"ax.hist(torque_segment_scores, bins=30, edgecolor='black', alpha=0.7, color='#9467bd')\n",
|
||
|
|
"ax.axvline(np.mean(torque_segment_scores), color='red', linestyle='--', linewidth=2, label=f'Mean = {np.mean(torque_segment_scores):.6f}')\n",
|
||
|
|
"ax.axvline(np.median(torque_segment_scores), color='orange', linestyle='--', linewidth=2, label=f'Median = {np.median(torque_segment_scores):.6f}')\n",
|
||
|
|
"ax.set_xlabel('R² Score', fontsize=12)\n",
|
||
|
|
"ax.set_ylabel('Number of Segments', fontsize=12)\n",
|
||
|
|
"ax.set_title('Distribution of Torque R² Scores Across All Parameter Segments', fontsize=14, fontweight='bold')\n",
|
||
|
|
"ax.legend(fontsize=11)\n",
|
||
|
|
"ax.grid(True, alpha=0.3, axis='y')\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "5bffc9a9",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Worst-Case Torque Prediction Visualization\n",
|
||
|
|
"\n",
|
||
|
|
"Let's visualize the worst-performing torque segment to assess model robustness."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "9d5580f3",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABKYAAAKyCAYAAADvidZRAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAA6AVJREFUeJzs/V1QXGmaJ3j+zwGBPkLgIIEUBMFITkTHxVh3VzqKvRjraattOZlX091WAREXa2M7uxWCLCvbvAirhJBNmpXFRbYSrDZ3q9KqM0GVa7O2e9ESVPVOrs1OZ+ARbTZbO/sh4VW1vX0RE4mHiiQJSUjgoNAHEpyzFyfdcQd3cMefc87zvv7/mYUFn6/+530euZyX9z3u+L7vg4iIiIiIiIiIKGJu3AGIiIiIiIiIiKg5cWGKiIiIiIiIiIhiwYUpIiIiIiIiIiKKBRemiIiIiIiIiIgoFlyYIiIiIiIiIiKiWHBhioiIiIiIiIiIYsGFKSIiIiIiIiIiigUXpoiIiIiIiIiIKBZcmCIiIlJgfn7+2N+by+UwPDyMrq4udHV1YXR0FLlcTjAdhS2TybBmRERE1JS4MEVETW98fByO41T9vOM4GBwcrPi5XC4Hx3GQyWTCincsmUymuEhRy382/kCcy+UwOjqKrq4uOI6Drq4uDA8PY3Z2Nu5oB4yPj+POnTvF97PZLBzHOfDf0NAQJicnkc/ny75/cHAQqVQKc3NzmJqaQjabxdDQ0IGvI726u7sxPDxsfc0ymQwcx1H597ARQ0NDcBwH2Ww27igNq/fxx0Ta6jU9PY2urq6q72uXzWYxPDyMoaEh6/5uE1E0uDBFRE1veHgYACouLhU+lsvlKj4ZL3w+nU6HF/AYrly5grm5ubL/3n//feTzeVy/fv3A55LJZNyRRc3Pz2NwcBC5XA43b97E4uIipqamAABzc3Mxpys3OzuLTCZTzFdqZmYGGxsb2NjYwNLSEq5fv45MJoPLly+XLSYWFqTS6TTGxsawsLCAfD6vbsEU2NvdVVjwrfeHmGw2W1xwHBwcxPj4eNnnM5lM8YfOwcFBTE9PVx2nNEelr6snaz6fLy6C7jc9PV3xB/3SBe9UKoWRkRGMjo7WOhV0TLlcDtPT0wd6Z//X1Fr7XC6HbDaLRCKBmZmZMCLHotbHn7AdVa96H1NsrddRjvvYW8v3Xbt2DQsLC1hcXMTc3JwVi5dEFDGfiKjJbWxs+AD8iYmJA58bGxvzR0ZGfAD+3Nzcgc+PjIz46XQ69Ixzc3P+xsZGQ2PMzMz4APzFxUWZUEotLi76APyRkZGKn290HqUlEgl/YWGh7GOFa6jUc77v+6lUyk+lUlXHXFhYUFnrwt+1kZERf2FhodiTU1NTNX3/3Nxc8e/q4uKiv7i4WDZHc3NzfiKR8GdmZvylpaXi+GNjY2XjFOZ3amrKX1pa8hcWFvxkMlnWM/VmHRsb85PJpF/pqdXU1JSfSCSKmQv/LS0tHfjaSv1gk0JvzszMxPZnJxKJIx8j6qn91NSUn06n/bGxMT+RSIR5CZFo9PFHSi31Os5jisZ6FR4jqr3fqOM+9tb6faX9kE6n1f07S0T6cWGKiMgPnlQlk8kDH08mk/7MzIyfSqUqPilOJBI1/1DdCAAN/7DaLAtTyWSyYi01mpmZqfjDx1E/GE5MTFRcAPH94Ie5VCpV9XvjNDY2duAH2qmpqarXUqrwA9JhCxrpdPrA35NKc5VKpaouVhUWi+rJWrrQddjCVC0q/bk2iXNhqtT+hchS9fZp4d+JwrWZvrDYyONPWKrV6ziPKRrrFfbC1HEfe2v9vsXFRT+VSvnpdFrlvz1EpB+P8hERITiKt/+4Xj6fRy6XQzqdRjqdPnAsKpvNIp/PqzvG18wKN5CenJyMO0pNpqam8P7779f9fZlM5sDxy8LRtFwuh8XFRYyMjEjFFJPJZPDBBx+UfWxsbKz4ucNMTk4imUwWv76ShYWFmv4+5nI5DA0NlX0slUoBQPGeM/VkvXbtGiYmJpBIJI78s48yOjqKbDZr5X3fTFFP7Qu1Kvw7Aeg7Liyt0uNPXOp9TGnGegHHf+yt9ftSqRQWFxexsLCg8t8eItKPC1NERKh8n6lMJoNEIoFkMlm8KXHpD4uFzxd+oAX2FgcK97/Zf9+a6enp4p81Pj6Orq6u4r2ACvfFKdyku3AvncI9awr3eDjsRu2NOG52YO+eToVXhJufn4fjOGXzNTk5eeBmroUbIe//IXx+fr44H0NDQzXfK6mwqFDrYs/+G6QPDQ0duBlutdo0mrXw5xfmtBaFnshms2X3RhkfH8fo6CjGx8eRTCaRyWSQyWRU3eej8Pdn/w+0iUQCiUTiyJsQ3759G+l0GvPz88U+reVm4fPz82V/R4HgHmyLi4tlHyv8+alUqq6ss7OzyOVyuH79+qE5gL2/N4fd+6rww3IUN2Uu/Tu/v7dr/fta7XHhsMcLzert01u3biGZTBa/Pp1Oi978WVONqj3+xOU4jylh1uuw+Tzq39cwHfext9HHbCKienBhiogIez8MLiwsFD9269at4gJH4fOliw77d2cUFjBSqRQ+++wzTE1NYWZm5sDNjAu7NTKZTPGH2eHhYXzwwQdYWlrC3NwcUqkU1tfXizfuBoLf6hZuRCvtuNkTiQTm5+cxOjqKZDKJubk5DA8P49q1a8fOMjs7W1xkWVxcxAcffFDcCXSUwivb1bpzZX5+Ht3d3Zibm8PS0hKSySSuXr1a/GEin89XrU2jWQu9tH/RpFRhYbL0VbHy+TwWFxeLvZfP54uLI6OjoxgeHi7+d/v27ZrmIQqFOatUm+7ubiwtLR36/YUF3Fu3bmFychI3b97E3bt3cfXq1QNfm8vlij0NAJ999lnZ52dmZnD79m1MTk4im80We3hqagrJZLLmrPl8HpOTk8W/C0flz+VymJubw+TkZPG/SpLJZNljURgKr9yYSCSKL4KQSqVw69atuseq9Lhw2Mc1q7dPZ2dny3aIFB4z5+fnG84Sd41qefyJ03EeU8KsF1B5Pmv99zUsx33sbfQxm4ioHq1xByAi0mL/cb1MJoObN28W30+lUlhYWCjbxr5/18rExETx1dVSqRRSqRQGBwcxPz9ffDJcOEZQWHAq/JkTExMAgh9KS5/0F35bWfgtZRiOmx0IjjEV5qZgaWnpWL8RzufzGB8fx8zMTHGeU6kUHj9+XHwyf9T316Mw5wU3b95EV1cXbt++jbGxMdy9e7fs60pr02jWwuLVYUdipqaminM/MzOD6elpzMzMlH1PIpGA7/v1XHaZ8fHxuo6NjY+PH+uoRiO7ZQr5Cj+gFyQSCQwPDyOTyZT9nSldHJyZmTnw9yaZTOKzzz7D1atXi306MTFRrHOtWScnJ9Hd3X2gj/ZLJBJIp9MV/45UWrBJJBLFHwrDMjo6ipGRkbL5PO5iQ6XHhcM+Xouo+nK/evq0cJy79KjT+++/j/Hxcdy6davhPHHXqJbHn4I46lXvY0rY9QIqz2et/74e5bhzfNzHXhN2OBKRPbgwRUT0W8PDw5icnEQ+n8f6+vqB+0eVbvnff2SscN+K/S9nnUwmi7/hLn3yWbpoceXKFQDA0NAQxsfHkU6nI71/RyPZC/flKjzhLhgeHj7WwlRhIWh8fPxAnsN2FpVmLuQ6zhwWFggKvwk+rDaNZq3lSX/pkZOpqSnMzs5icnJS9J4oGo7kALXtctt/r5PSY2+lf1cL9SscPZqbmytbFCrskJqbm8PIyEjZbrNadioVjrHMzs7W9PVjY2MH7o1V+Dty9+7dA4sN3d3dod5jKpfLIZfLida+2ljH/TO09OV+pX1aWPRMJpNlf59TqVTDO3A01Kiexx+N9dr/mBJmvfb/OQX1/vta67iSjvsLLxN2QBKROXiUj4jotwpPDgv350kmk2VPvD744APk83lks9kDny8sUnR3dx8YN5lMHvghc/+Ol8JvV8fHxzE4OFjTvXOkNJK9ll0/x7GxsQE/eOXY4n+17LqodK+wo2QyGYyOjhbvkVWqltocN+txTE1NYX5+XvWNsQv32yr9L5vNFnukUl8
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1200x700 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Worst-case torque prediction achieves R² = -0.638045\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Following generated by AI - see [11]\n",
|
||
|
|
"# Get worst torque segment data\n",
|
||
|
|
"worst_torque_data = magDf.iloc[worst_torque_segment['start']:worst_torque_segment['end']]\n",
|
||
|
|
"worst_torque_currL = worst_torque_segment['currL']\n",
|
||
|
|
"worst_torque_currR = worst_torque_segment['currR']\n",
|
||
|
|
"worst_torque_roll = worst_torque_segment['roll']\n",
|
||
|
|
"worst_torque_r2 = worst_torque_segment['r2']\n",
|
||
|
|
"\n",
|
||
|
|
"# Create fine grid for this worst segment\n",
|
||
|
|
"gap_fine_worst_torque = np.linspace(2, 30, 500)\n",
|
||
|
|
"X_worst_torque_fine = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [worst_torque_currL] * 500,\n",
|
||
|
|
" 'currR [A]': [worst_torque_currR] * 500,\n",
|
||
|
|
" 'rollDeg [deg]': [worst_torque_roll] * 500,\n",
|
||
|
|
" 'invGap': 1/gap_fine_worst_torque\n",
|
||
|
|
"})\n",
|
||
|
|
"\n",
|
||
|
|
"# Get predictions\n",
|
||
|
|
"X_worst_torque_poly = poly_best.transform(X_worst_torque_fine)\n",
|
||
|
|
"y_worst_torque_pred = model_best.predict(X_worst_torque_poly)\n",
|
||
|
|
"\n",
|
||
|
|
"# Extract actual data\n",
|
||
|
|
"gap_worst_torque_actual = worst_torque_data['GapHeight [mm]'].values\n",
|
||
|
|
"torque_worst_actual = worst_torque_data['YokeTorque.Torque [mNewtonMeter]'].values\n",
|
||
|
|
"\n",
|
||
|
|
"# Create the plot\n",
|
||
|
|
"fig, ax = plt.subplots(1, 1, figsize=(12, 7))\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot polynomial curve for TORQUE\n",
|
||
|
|
"ax.plot(gap_fine_worst_torque, y_worst_torque_pred[:, 1], '-', linewidth=2.5, \n",
|
||
|
|
" label=f'Degree {best_degree} Polynomial', color='#ff7f0e', alpha=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot actual data points\n",
|
||
|
|
"ax.scatter(gap_worst_torque_actual, torque_worst_actual, s=120, marker='o', \n",
|
||
|
|
" color='#9467bd', edgecolors='black', linewidths=1.5,\n",
|
||
|
|
" label='Ansys Data', zorder=5)\n",
|
||
|
|
"\n",
|
||
|
|
"# Formatting\n",
|
||
|
|
"ax.set_xlabel('Gap Height (mm)', fontsize=13)\n",
|
||
|
|
"ax.set_ylabel('Torque (mN·m)', fontsize=13)\n",
|
||
|
|
"ax.set_title(f'Worst Torque Case (R² = {worst_torque_r2:.6f}) - currL={worst_torque_currL:.0f}A, currR={worst_torque_currR:.0f}A, roll={worst_torque_roll:.1f}°', \n",
|
||
|
|
" fontsize=14, fontweight='bold')\n",
|
||
|
|
"ax.legend(fontsize=12, loc='best')\n",
|
||
|
|
"ax.grid(True, alpha=0.3)\n",
|
||
|
|
"\n",
|
||
|
|
"# Add vertical lines at data points\n",
|
||
|
|
"for gap_point in gap_worst_torque_actual:\n",
|
||
|
|
" ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.3, linewidth=0.8)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Worst-case torque prediction achieves R² = {worst_torque_r2:.6f}\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "c63aca1e",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### We can see here that at roll angle and current differential of 0 (when torque will be very close to 0):\n",
|
||
|
|
"\n",
|
||
|
|
"There is significant overfitting to noise, but within the operating range of 6 - 24 mm, there is never more variation than there would be with simple linear interpolation.\n",
|
||
|
|
"\n",
|
||
|
|
"Hopefully this means we are ok"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "1e1c384f",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Export Model Summary\n",
|
||
|
|
"\n",
|
||
|
|
"Display key model information before exporting."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "778f58df",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"============================================================\n",
|
||
|
|
"MODEL EXPORT SUMMARY\n",
|
||
|
|
"============================================================\n",
|
||
|
|
"\n",
|
||
|
|
"Model Type: Polynomial Regression\n",
|
||
|
|
"Polynomial Degree: 6\n",
|
||
|
|
"Number of Features: 4\n",
|
||
|
|
"Number of Polynomial Terms: 210\n",
|
||
|
|
"\n",
|
||
|
|
"Input Features:\n",
|
||
|
|
" - currL [A]\n",
|
||
|
|
" - currR [A]\n",
|
||
|
|
" - rollDeg [deg]\n",
|
||
|
|
" - invGap (1/GapHeight)\n",
|
||
|
|
"\n",
|
||
|
|
"Outputs:\n",
|
||
|
|
" - YokeForce.Force_z [newton]\n",
|
||
|
|
" - YokeTorque.Torque [mNewtonMeter]\n",
|
||
|
|
"\n",
|
||
|
|
"Model Performance:\n",
|
||
|
|
" Force R² Score: 0.999967\n",
|
||
|
|
" Torque R² Score: 0.999895\n",
|
||
|
|
"\n",
|
||
|
|
"Model Coefficients Shape:\n",
|
||
|
|
" Force coefficients: (210,)\n",
|
||
|
|
" Torque coefficients: (210,)\n",
|
||
|
|
"============================================================\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was generated by AI - see [12]\n",
|
||
|
|
"print(\"=\" * 60)\n",
|
||
|
|
"print(\"MODEL EXPORT SUMMARY\")\n",
|
||
|
|
"print(\"=\" * 60)\n",
|
||
|
|
"print(f\"\\nModel Type: Polynomial Regression\")\n",
|
||
|
|
"print(f\"Polynomial Degree: {best_degree}\")\n",
|
||
|
|
"print(f\"Number of Features: {poly_best.n_features_in_}\")\n",
|
||
|
|
"print(f\"Number of Polynomial Terms: {poly_best.n_output_features_}\")\n",
|
||
|
|
"print(f\"\\nInput Features:\")\n",
|
||
|
|
"print(\" - currL [A]\")\n",
|
||
|
|
"print(\" - currR [A]\")\n",
|
||
|
|
"print(\" - rollDeg [deg]\")\n",
|
||
|
|
"print(\" - invGap (1/GapHeight)\")\n",
|
||
|
|
"print(f\"\\nOutputs:\")\n",
|
||
|
|
"print(\" - YokeForce.Force_z [newton]\")\n",
|
||
|
|
"print(\" - YokeTorque.Torque [mNewtonMeter]\")\n",
|
||
|
|
"print(f\"\\nModel Performance:\")\n",
|
||
|
|
"print(f\" Force R² Score: {r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]):.6f}\")\n",
|
||
|
|
"print(f\" Torque R² Score: {torque_r2_overall:.6f}\")\n",
|
||
|
|
"print(f\"\\nModel Coefficients Shape:\")\n",
|
||
|
|
"print(f\" Force coefficients: {model_best.coef_[0].shape}\")\n",
|
||
|
|
"print(f\" Torque coefficients: {model_best.coef_[1].shape}\")\n",
|
||
|
|
"print(\"=\" * 60)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "88208b13",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Create Standalone Python Predictor Class\n",
|
||
|
|
"\n",
|
||
|
|
"Generate a self-contained Python module that can be imported into your simulator without scikit-learn dependencies."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "edb239e1",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"✓ Saved standalone predictor: maglev_predictor.py\n",
|
||
|
|
"\n",
|
||
|
|
"This module:\n",
|
||
|
|
" - Has NO scikit-learn dependency (only numpy)\n",
|
||
|
|
" - Can be imported directly into your simulator\n",
|
||
|
|
" - Includes both single and batch prediction methods\n",
|
||
|
|
" - Implements polynomial feature generation internally\n",
|
||
|
|
"\n",
|
||
|
|
"To use:\n",
|
||
|
|
"```python\n",
|
||
|
|
"from maglev_predictor import MaglevPredictor\n",
|
||
|
|
"predictor = MaglevPredictor()\n",
|
||
|
|
"force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)\n",
|
||
|
|
"```\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was generated by AI - see [12]\n",
|
||
|
|
"# Generate standalone predictor class\n",
|
||
|
|
"predictor_code = f'''\"\"\"\n",
|
||
|
|
"Magnetic Levitation Force and Torque Predictor\n",
|
||
|
|
"Generated from polynomial regression model (degree {best_degree})\n",
|
||
|
|
"\n",
|
||
|
|
"Performance:\n",
|
||
|
|
" - Force R²: {r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]):.6f}\n",
|
||
|
|
" - Torque R²: {torque_r2_overall:.6f}\n",
|
||
|
|
"\n",
|
||
|
|
"Usage:\n",
|
||
|
|
" predictor = MaglevPredictor()\n",
|
||
|
|
" force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"from itertools import combinations_with_replacement\n",
|
||
|
|
"\n",
|
||
|
|
"class MaglevPredictor:\n",
|
||
|
|
" def __init__(self):\n",
|
||
|
|
" \"\"\"Initialize the magnetic levitation force/torque predictor.\"\"\"\n",
|
||
|
|
" self.degree = {best_degree}\n",
|
||
|
|
" self.n_features = 4 # currL, currR, roll, invGap\n",
|
||
|
|
" \n",
|
||
|
|
" # Force model coefficients\n",
|
||
|
|
" self.force_intercept = {model_best.intercept_[0]}\n",
|
||
|
|
" self.force_coef = np.array({model_best.coef_[0].tolist()})\n",
|
||
|
|
" \n",
|
||
|
|
" # Torque model coefficients \n",
|
||
|
|
" self.torque_intercept = {model_best.intercept_[1]}\n",
|
||
|
|
" self.torque_coef = np.array({model_best.coef_[1].tolist()})\n",
|
||
|
|
" \n",
|
||
|
|
" def _polynomial_features(self, X):\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" Generate polynomial features up to specified degree.\n",
|
||
|
|
" Mimics sklearn's PolynomialFeatures with include_bias=True.\n",
|
||
|
|
" \n",
|
||
|
|
" Args:\n",
|
||
|
|
" X: numpy array of shape (n_samples, 4) with [currL, currR, roll, invGap]\n",
|
||
|
|
" \n",
|
||
|
|
" Returns:\n",
|
||
|
|
" Polynomial features array\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" n_samples = X.shape[0]\n",
|
||
|
|
" \n",
|
||
|
|
" # Start with bias term (column of ones)\n",
|
||
|
|
" features = [np.ones(n_samples)]\n",
|
||
|
|
" \n",
|
||
|
|
" # Add original features\n",
|
||
|
|
" for i in range(self.n_features):\n",
|
||
|
|
" features.append(X[:, i])\n",
|
||
|
|
" \n",
|
||
|
|
" # Add polynomial combinations\n",
|
||
|
|
" for deg in range(2, self.degree + 1):\n",
|
||
|
|
" for combo in combinations_with_replacement(range(self.n_features), deg):\n",
|
||
|
|
" term = np.ones(n_samples)\n",
|
||
|
|
" for idx in combo:\n",
|
||
|
|
" term *= X[:, idx]\n",
|
||
|
|
" features.append(term)\n",
|
||
|
|
" \n",
|
||
|
|
" return np.column_stack(features)\n",
|
||
|
|
" \n",
|
||
|
|
" def predict(self, currL, currR, roll, gap_height):\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" Predict force and torque for given operating conditions.\n",
|
||
|
|
" \n",
|
||
|
|
" Args:\n",
|
||
|
|
" currL: Left coil current in Amps\n",
|
||
|
|
" currR: Right coil current in Amps\n",
|
||
|
|
" roll: Roll angle in degrees\n",
|
||
|
|
" gap_height: Gap height in mm\n",
|
||
|
|
" \n",
|
||
|
|
" Returns:\n",
|
||
|
|
" tuple: (force [N], torque [mN·m])\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" # Compute inverse gap (critical transformation!)\n",
|
||
|
|
" invGap = 1.0 / gap_height\n",
|
||
|
|
" \n",
|
||
|
|
" # Create input array\n",
|
||
|
|
" X = np.array([[currL, currR, roll, invGap]])\n",
|
||
|
|
" \n",
|
||
|
|
" # Generate polynomial features\n",
|
||
|
|
" X_poly = self._polynomial_features(X)\n",
|
||
|
|
" \n",
|
||
|
|
" # Compute predictions\n",
|
||
|
|
" force = self.force_intercept + np.dot(X_poly, self.force_coef)[0]\n",
|
||
|
|
" torque = self.torque_intercept + np.dot(X_poly, self.torque_coef)[0]\n",
|
||
|
|
" \n",
|
||
|
|
" return force, torque\n",
|
||
|
|
" \n",
|
||
|
|
" def predict_batch(self, currL_array, currR_array, roll_array, gap_height_array):\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" Predict force and torque for multiple operating conditions.\n",
|
||
|
|
" \n",
|
||
|
|
" Args:\n",
|
||
|
|
" currL_array: Array of left coil currents [A]\n",
|
||
|
|
" currR_array: Array of right coil currents [A]\n",
|
||
|
|
" roll_array: Array of roll angles [deg]\n",
|
||
|
|
" gap_height_array: Array of gap heights [mm]\n",
|
||
|
|
" \n",
|
||
|
|
" Returns:\n",
|
||
|
|
" tuple: (force_array [N], torque_array [mN·m])\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" # Convert to numpy arrays\n",
|
||
|
|
" currL_array = np.asarray(currL_array)\n",
|
||
|
|
" currR_array = np.asarray(currR_array)\n",
|
||
|
|
" roll_array = np.asarray(roll_array)\n",
|
||
|
|
" gap_height_array = np.asarray(gap_height_array)\n",
|
||
|
|
" \n",
|
||
|
|
" # Compute inverse gaps\n",
|
||
|
|
" invGap_array = 1.0 / gap_height_array\n",
|
||
|
|
" \n",
|
||
|
|
" # Stack into feature matrix\n",
|
||
|
|
" X = np.column_stack([currL_array, currR_array, roll_array, invGap_array])\n",
|
||
|
|
" \n",
|
||
|
|
" # Generate polynomial features\n",
|
||
|
|
" X_poly = self._polynomial_features(X)\n",
|
||
|
|
" \n",
|
||
|
|
" # Compute predictions\n",
|
||
|
|
" force_array = self.force_intercept + np.dot(X_poly, self.force_coef)\n",
|
||
|
|
" torque_array = self.torque_intercept + np.dot(X_poly, self.torque_coef)\n",
|
||
|
|
" \n",
|
||
|
|
" return force_array, torque_array\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" # Example usage\n",
|
||
|
|
" predictor = MaglevPredictor()\n",
|
||
|
|
" \n",
|
||
|
|
" # Single prediction\n",
|
||
|
|
" force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)\n",
|
||
|
|
" print(f\"Single prediction:\")\n",
|
||
|
|
" print(f\" Force: {{force:.2f}} N\")\n",
|
||
|
|
" print(f\" Torque: {{torque:.2f}} mN·m\")\n",
|
||
|
|
" \n",
|
||
|
|
" # Batch prediction\n",
|
||
|
|
" currL = np.array([-15, -15, -10])\n",
|
||
|
|
" currR = np.array([-15, -10, -10])\n",
|
||
|
|
" roll = np.array([0, 0.5, 1.0])\n",
|
||
|
|
" gap = np.array([10, 12, 15])\n",
|
||
|
|
" \n",
|
||
|
|
" forces, torques = predictor.predict_batch(currL, currR, roll, gap)\n",
|
||
|
|
" print(f\"\\\\nBatch prediction:\")\n",
|
||
|
|
" for i in range(len(forces)):\n",
|
||
|
|
" print(f\" Condition {{i+1}}: Force={{forces[i]:.2f}} N, Torque={{torques[i]:.2f}} mN·m\")\n",
|
||
|
|
"'''\n",
|
||
|
|
"\n",
|
||
|
|
"# Save to file\n",
|
||
|
|
"predictor_filename = 'maglev_predictor.py'\n",
|
||
|
|
"with open(predictor_filename, 'w') as f:\n",
|
||
|
|
" f.write(predictor_code)\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"✓ Saved standalone predictor: {predictor_filename}\")\n",
|
||
|
|
"print(f\"\\nThis module:\")\n",
|
||
|
|
"print(f\" - Has NO scikit-learn dependency (only numpy)\")\n",
|
||
|
|
"print(f\" - Can be imported directly into your simulator\")\n",
|
||
|
|
"print(f\" - Includes both single and batch prediction methods\")\n",
|
||
|
|
"print(f\" - Implements polynomial feature generation internally\")\n",
|
||
|
|
"print(f\"\\nTo use:\")\n",
|
||
|
|
"print(f\"```python\")\n",
|
||
|
|
"print(f\"from maglev_predictor import MaglevPredictor\")\n",
|
||
|
|
"print(f\"predictor = MaglevPredictor()\")\n",
|
||
|
|
"print(f\"force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)\")\n",
|
||
|
|
"print(f\"```\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "e4f40fb4",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Test the Exported Predictor\n",
|
||
|
|
"\n",
|
||
|
|
"Verify the standalone predictor produces identical results to the original model."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "4a1ddacd",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Validation: Comparing Standalone Predictor vs Original Model\n",
|
||
|
|
"================================================================================\n",
|
||
|
|
"\n",
|
||
|
|
"Test Case 1: currL=-15A, currR=-15A, roll=0.0°, gap=10.0mm\n",
|
||
|
|
" Force: Standalone= 103.3304 N | Original= 103.3304 N | Diff=2.84e-14\n",
|
||
|
|
" Torque: Standalone= -0.8561 mN·m | Original= -0.8561 mN·m | Diff=1.02e-12\n",
|
||
|
|
"\n",
|
||
|
|
"Test Case 2: currL=-15A, currR=-15A, roll=1.0°, gap=10.0mm\n",
|
||
|
|
" Force: Standalone= 104.0520 N | Original= 104.0520 N | Diff=2.84e-14\n",
|
||
|
|
" Torque: Standalone= 450.3916 mN·m | Original= 450.3916 mN·m | Diff=3.41e-13\n",
|
||
|
|
"\n",
|
||
|
|
"Test Case 3: currL=-15A, currR=-10A, roll=0.5°, gap=12.0mm\n",
|
||
|
|
" Force: Standalone= 75.0401 N | Original= 75.0401 N | Diff=5.68e-14\n",
|
||
|
|
" Torque: Standalone= 390.7023 mN·m | Original= 390.7023 mN·m | Diff=3.98e-13\n",
|
||
|
|
"\n",
|
||
|
|
"Test Case 4: currL=-10A, currR=-10A, roll=2.0°, gap=15.0mm\n",
|
||
|
|
" Force: Standalone= 50.4344 N | Original= 50.4344 N | Diff=1.42e-14\n",
|
||
|
|
" Torque: Standalone= 303.7400 mN·m | Original= 303.7400 mN·m | Diff=1.71e-13\n",
|
||
|
|
"\n",
|
||
|
|
"================================================================================\n",
|
||
|
|
"✓ All tests passed! Standalone predictor matches original model perfectly.\n",
|
||
|
|
"\n",
|
||
|
|
"The standalone predictor is ready for integration into your simulator!\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was generated by AI - see [12]\n",
|
||
|
|
"# Import the standalone predictor (reload to get latest version)\n",
|
||
|
|
"import importlib\n",
|
||
|
|
"import maglev_predictor\n",
|
||
|
|
"importlib.reload(maglev_predictor)\n",
|
||
|
|
"from maglev_predictor import MaglevPredictor\n",
|
||
|
|
"\n",
|
||
|
|
"# Create predictor instance\n",
|
||
|
|
"standalone_predictor = MaglevPredictor()\n",
|
||
|
|
"\n",
|
||
|
|
"# Test cases\n",
|
||
|
|
"test_cases = [\n",
|
||
|
|
" {'currL': -15, 'currR': -15, 'roll': 0.0, 'gap': 10.0},\n",
|
||
|
|
" {'currL': -15, 'currR': -15, 'roll': 1.0, 'gap': 10.0},\n",
|
||
|
|
" {'currL': -15, 'currR': -10, 'roll': 0.5, 'gap': 12.0},\n",
|
||
|
|
" {'currL': -10, 'currR': -10, 'roll': 2.0, 'gap': 15.0},\n",
|
||
|
|
"]\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"Validation: Comparing Standalone Predictor vs Original Model\")\n",
|
||
|
|
"print(\"=\" * 80)\n",
|
||
|
|
"\n",
|
||
|
|
"for i, case in enumerate(test_cases):\n",
|
||
|
|
" # Standalone predictor\n",
|
||
|
|
" force_standalone, torque_standalone = standalone_predictor.predict(\n",
|
||
|
|
" case['currL'], case['currR'], case['roll'], case['gap']\n",
|
||
|
|
" )\n",
|
||
|
|
" \n",
|
||
|
|
" # Original model\n",
|
||
|
|
" X_test = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': [case['currL']],\n",
|
||
|
|
" 'currR [A]': [case['currR']],\n",
|
||
|
|
" 'rollDeg [deg]': [case['roll']],\n",
|
||
|
|
" 'invGap': [1/case['gap']]\n",
|
||
|
|
" })\n",
|
||
|
|
" X_test_poly = poly_best.transform(X_test)\n",
|
||
|
|
" y_test_pred = model_best.predict(X_test_poly)\n",
|
||
|
|
" force_original = y_test_pred[0, 0]\n",
|
||
|
|
" torque_original = y_test_pred[0, 1]\n",
|
||
|
|
" \n",
|
||
|
|
" # Compare\n",
|
||
|
|
" force_diff = abs(force_standalone - force_original)\n",
|
||
|
|
" torque_diff = abs(torque_standalone - torque_original)\n",
|
||
|
|
" \n",
|
||
|
|
" print(f\"\\nTest Case {i+1}: currL={case['currL']}A, currR={case['currR']}A, \" + \n",
|
||
|
|
" f\"roll={case['roll']}°, gap={case['gap']}mm\")\n",
|
||
|
|
" print(f\" Force: Standalone={force_standalone:9.4f} N | Original={force_original:9.4f} N | Diff={force_diff:.2e}\")\n",
|
||
|
|
" print(f\" Torque: Standalone={torque_standalone:9.4f} mN·m | Original={torque_original:9.4f} mN·m | Diff={torque_diff:.2e}\")\n",
|
||
|
|
" \n",
|
||
|
|
" # Verify match (should be essentially identical, accounting for floating point)\n",
|
||
|
|
" assert force_diff < 1e-8, f\"Force mismatch in test case {i+1}\"\n",
|
||
|
|
" assert torque_diff < 1e-8, f\"Torque mismatch in test case {i+1}\"\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"\\n\" + \"=\" * 80)\n",
|
||
|
|
"print(\"✓ All tests passed! Standalone predictor matches original model perfectly.\")\n",
|
||
|
|
"print(\"\\nThe standalone predictor is ready for integration into your simulator!\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "d473d5fa",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Export Sklearn Model (Lossless & Faster)\n",
|
||
|
|
"\n",
|
||
|
|
"Save the actual sklearn model using joblib for lossless serialization. This is also faster than the standalone predictor"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "5fd33470",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"✓ Saved sklearn model: maglev_model.pkl\n",
|
||
|
|
"\n",
|
||
|
|
"File size: 6.13 KB\n",
|
||
|
|
"\n",
|
||
|
|
"Contents:\n",
|
||
|
|
" - PolynomialFeatures transformer (degree 6)\n",
|
||
|
|
" - LinearRegression model\n",
|
||
|
|
" - Feature names and metadata\n",
|
||
|
|
" - Performance metrics\n",
|
||
|
|
"\n",
|
||
|
|
"To use:\n",
|
||
|
|
"```python\n",
|
||
|
|
"import joblib\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import pandas as pd\n",
|
||
|
|
"\n",
|
||
|
|
"# Load model\n",
|
||
|
|
"data = joblib.load('maglev_model.pkl')\n",
|
||
|
|
"poly = data['poly_features']\n",
|
||
|
|
"model = data['model']\n",
|
||
|
|
"\n",
|
||
|
|
"# Make prediction\n",
|
||
|
|
"X = pd.DataFrame({'currL [A]': [-15], 'currR [A]': [-15], \n",
|
||
|
|
" 'rollDeg [deg]': [1.0], 'invGap': [1/10.0]})\n",
|
||
|
|
"X_poly = poly.transform(X)\n",
|
||
|
|
"force, torque = model.predict(X_poly)[0]\n",
|
||
|
|
"```\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was generated by AI - see [12]\n",
|
||
|
|
"import joblib\n",
|
||
|
|
"import pickle\n",
|
||
|
|
"\n",
|
||
|
|
"# Save the polynomial transformer and model\n",
|
||
|
|
"model_data = {\n",
|
||
|
|
" 'poly_features': poly_best,\n",
|
||
|
|
" 'model': model_best,\n",
|
||
|
|
" 'degree': best_degree,\n",
|
||
|
|
" 'feature_names': ['currL [A]', 'currR [A]', 'rollDeg [deg]', 'invGap'],\n",
|
||
|
|
" 'performance': {\n",
|
||
|
|
" 'force_r2': r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]),\n",
|
||
|
|
" 'torque_r2': torque_r2_overall\n",
|
||
|
|
" }\n",
|
||
|
|
"}\n",
|
||
|
|
"\n",
|
||
|
|
"# Save using joblib (preferred for sklearn models)\n",
|
||
|
|
"joblib.dump(model_data, 'maglev_model.pkl')\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"✓ Saved sklearn model: maglev_model.pkl\")\n",
|
||
|
|
"print(f\"\\nFile size: {os.path.getsize('maglev_model.pkl') / 1024:.2f} KB\")\n",
|
||
|
|
"print(f\"\\nContents:\")\n",
|
||
|
|
"print(f\" - PolynomialFeatures transformer (degree {best_degree})\")\n",
|
||
|
|
"print(f\" - LinearRegression model\")\n",
|
||
|
|
"print(f\" - Feature names and metadata\")\n",
|
||
|
|
"print(f\" - Performance metrics\")\n",
|
||
|
|
"print(f\"\\nTo use:\")\n",
|
||
|
|
"print(f\"```python\")\n",
|
||
|
|
"print(f\"import joblib\")\n",
|
||
|
|
"print(f\"import numpy as np\")\n",
|
||
|
|
"print(f\"import pandas as pd\")\n",
|
||
|
|
"print(f\"\")\n",
|
||
|
|
"print(f\"# Load model\")\n",
|
||
|
|
"print(f\"data = joblib.load('maglev_model.pkl')\")\n",
|
||
|
|
"print(f\"poly = data['poly_features']\")\n",
|
||
|
|
"print(f\"model = data['model']\")\n",
|
||
|
|
"print(f\"\")\n",
|
||
|
|
"print(f\"# Make prediction\")\n",
|
||
|
|
"print(f\"X = pd.DataFrame({{'currL [A]': [-15], 'currR [A]': [-15], \")\n",
|
||
|
|
"print(f\" 'rollDeg [deg]': [1.0], 'invGap': [1/10.0]}})\")\n",
|
||
|
|
"print(f\"X_poly = poly.transform(X)\")\n",
|
||
|
|
"print(f\"force, torque = model.predict(X_poly)[0]\")\n",
|
||
|
|
"print(f\"```\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "99b93e74",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Compare Performance: Sklearn Model vs Standalone Predictor\n",
|
||
|
|
"\n",
|
||
|
|
"Benchmark both approaches to show the speed advantage of using sklearn directly."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "889df820",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Performance Comparison (1000 predictions)\n",
|
||
|
|
"============================================================\n",
|
||
|
|
"Sklearn Model (pickled): 2.48 ms\n",
|
||
|
|
"Standalone Predictor: 4.07 ms\n",
|
||
|
|
"\n",
|
||
|
|
"Speedup: 1.64x faster with sklearn model\n",
|
||
|
|
"\n",
|
||
|
|
"Result verification:\n",
|
||
|
|
" Max force difference: 2.27e-13 N\n",
|
||
|
|
" Max torque difference: 5.46e-12 mN·m\n",
|
||
|
|
" ✓ Results match!\n",
|
||
|
|
"\n",
|
||
|
|
"Per-prediction time:\n",
|
||
|
|
" Sklearn: 2.48 μs\n",
|
||
|
|
" Standalone: 4.07 μs\n",
|
||
|
|
"\n",
|
||
|
|
"============================================================\n",
|
||
|
|
"Recommendation:\n",
|
||
|
|
" Use sklearn model (maglev_model.pkl) - 1.6x faster!\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was generated by AI - see [12]\n",
|
||
|
|
"import time\n",
|
||
|
|
"\n",
|
||
|
|
"# Load the sklearn model\n",
|
||
|
|
"loaded_data = joblib.load('maglev_model.pkl')\n",
|
||
|
|
"loaded_poly = loaded_data['poly_features']\n",
|
||
|
|
"loaded_model = loaded_data['model']\n",
|
||
|
|
"\n",
|
||
|
|
"# Create test data (1000 predictions)\n",
|
||
|
|
"n_tests = 1000\n",
|
||
|
|
"test_currL = np.random.uniform(-15, 0, n_tests)\n",
|
||
|
|
"test_currR = np.random.uniform(-15, 0, n_tests)\n",
|
||
|
|
"test_roll = np.random.uniform(-5, 5, n_tests)\n",
|
||
|
|
"test_gap = np.random.uniform(6, 24, n_tests)\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"Performance Comparison (1000 predictions)\")\n",
|
||
|
|
"print(\"=\" * 60)\n",
|
||
|
|
"\n",
|
||
|
|
"# Benchmark 1: Sklearn model (loaded from pickle)\n",
|
||
|
|
"start = time.perf_counter()\n",
|
||
|
|
"X_test = pd.DataFrame({\n",
|
||
|
|
" 'currL [A]': test_currL,\n",
|
||
|
|
" 'currR [A]': test_currR,\n",
|
||
|
|
" 'rollDeg [deg]': test_roll,\n",
|
||
|
|
" 'invGap': 1/test_gap\n",
|
||
|
|
"})\n",
|
||
|
|
"X_test_poly = loaded_poly.transform(X_test)\n",
|
||
|
|
"sklearn_results = loaded_model.predict(X_test_poly)\n",
|
||
|
|
"sklearn_time = time.perf_counter() - start\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Sklearn Model (pickled): {sklearn_time*1000:.2f} ms\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Benchmark 2: Standalone predictor\n",
|
||
|
|
"start = time.perf_counter()\n",
|
||
|
|
"standalone_results = standalone_predictor.predict_batch(test_currL, test_currR, test_roll, test_gap)\n",
|
||
|
|
"standalone_time = time.perf_counter() - start\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"Standalone Predictor: {standalone_time*1000:.2f} ms\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Speedup\n",
|
||
|
|
"speedup = standalone_time / sklearn_time\n",
|
||
|
|
"print(f\"\\nSpeedup: {speedup:.2f}x {'faster' if speedup > 1 else 'slower'} with sklearn model\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Verify results match\n",
|
||
|
|
"force_diff = np.abs(sklearn_results[:, 0] - standalone_results[0]).max()\n",
|
||
|
|
"torque_diff = np.abs(sklearn_results[:, 1] - standalone_results[1]).max()\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"\\nResult verification:\")\n",
|
||
|
|
"print(f\" Max force difference: {force_diff:.2e} N\")\n",
|
||
|
|
"print(f\" Max torque difference: {torque_diff:.2e} mN·m\")\n",
|
||
|
|
"print(f\" ✓ Results match!\" if (force_diff < 1e-8 and torque_diff < 1e-8) else \" ✗ Results differ!\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Per-prediction timing\n",
|
||
|
|
"print(f\"\\nPer-prediction time:\")\n",
|
||
|
|
"print(f\" Sklearn: {sklearn_time/n_tests*1e6:.2f} μs\")\n",
|
||
|
|
"print(f\" Standalone: {standalone_time/n_tests*1e6:.2f} μs\")\n",
|
||
|
|
"\n",
|
||
|
|
"print(f\"\\n{'='*60}\")\n",
|
||
|
|
"print(f\"Recommendation:\")\n",
|
||
|
|
"if speedup > 1.2:\n",
|
||
|
|
" print(f\" Use sklearn model (maglev_model.pkl) - {speedup:.1f}x faster!\")\n",
|
||
|
|
"elif speedup < 0.8:\n",
|
||
|
|
" print(f\" Use standalone predictor - {1/speedup:.1f}x faster!\")\n",
|
||
|
|
"else:\n",
|
||
|
|
" print(f\" Performance is similar - choose based on dependencies:\")\n",
|
||
|
|
" print(f\" • Sklearn model: Faster, requires sklearn+joblib\")\n",
|
||
|
|
" print(f\" • Standalone: No sklearn dependency, portable\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "21babeb3",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"### Now, let's find the equilibrium height for our pod, given mass of 5.8 kg. \n",
|
||
|
|
"\n",
|
||
|
|
"5.8 kg * 9.81 $m/s^2$ = 56.898 N"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "badbc379",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"======================================================================\n",
|
||
|
|
"EQUILIBRIUM GAP HEIGHT FINDER (Analytical Solution)\n",
|
||
|
|
"======================================================================\n",
|
||
|
|
"Pod mass: 5.8 kg\n",
|
||
|
|
"Total weight: 56.898 N\n",
|
||
|
|
"Target force per yoke: 28.449 N\n",
|
||
|
|
"Parameters: currL = 0 A, currR = 0 A, roll = 0°\n",
|
||
|
|
"\n",
|
||
|
|
"using scipy.optimize.fsolve:\n",
|
||
|
|
" Gap: 16.491741 mm → Force: 28.449000 N\n",
|
||
|
|
"\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The following was generated by AI - see [13]\n",
|
||
|
|
"# Find equilibrium gap height for 5.8 kg pod using polynomial root finding\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"from maglev_predictor import MaglevPredictor\n",
|
||
|
|
"from scipy.optimize import fsolve\n",
|
||
|
|
"\n",
|
||
|
|
"# Initialize predictor\n",
|
||
|
|
"predictor = MaglevPredictor()\n",
|
||
|
|
"\n",
|
||
|
|
"# Target force for 5.8 kg pod (total force = weight)\n",
|
||
|
|
"# Since we have TWO yokes (front and back), each produces this force\n",
|
||
|
|
"target_force_per_yoke = 5.8 * 9.81 / 2 # 28.449 N per yoke\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"=\" * 70)\n",
|
||
|
|
"print(\"EQUILIBRIUM GAP HEIGHT FINDER (Analytical Solution)\")\n",
|
||
|
|
"print(\"=\" * 70)\n",
|
||
|
|
"print(f\"Pod mass: 5.8 kg\")\n",
|
||
|
|
"print(f\"Total weight: {5.8 * 9.81:.3f} N\")\n",
|
||
|
|
"print(f\"Target force per yoke: {target_force_per_yoke:.3f} N\")\n",
|
||
|
|
"print(f\"Parameters: currL = 0 A, currR = 0 A, roll = 0°\")\n",
|
||
|
|
"print()\n",
|
||
|
|
"\n",
|
||
|
|
"# Method 2: Use scipy.optimize.fsolve for verification\n",
|
||
|
|
"def force_error(gap_height):\n",
|
||
|
|
" # Handle array input from fsolve (convert to scalar)\n",
|
||
|
|
" gap_height = float(np.atleast_1d(gap_height)[0])\n",
|
||
|
|
" force, _ = predictor.predict(currL=0, currR=0, roll=0, gap_height=gap_height)\n",
|
||
|
|
" return force - target_force_per_yoke\n",
|
||
|
|
"\n",
|
||
|
|
"# Try multiple initial guesses to find all solutions\n",
|
||
|
|
"initial_guesses = [8, 10, 15, 20, 25]\n",
|
||
|
|
"scipy_solutions = []\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"using scipy.optimize.fsolve:\")\n",
|
||
|
|
"for guess in initial_guesses:\n",
|
||
|
|
" solution = fsolve(force_error, guess)[0]\n",
|
||
|
|
" if 5 <= solution <= 30: # Valid range\n",
|
||
|
|
" force, torque = predictor.predict(currL=0, currR=0, roll=0, gap_height=solution)\n",
|
||
|
|
" error = abs(force - target_force_per_yoke)\n",
|
||
|
|
" if error < 0.01: # Valid solution (within 10 mN)\n",
|
||
|
|
" scipy_solutions.append((solution, force))\n",
|
||
|
|
"\n",
|
||
|
|
"# Remove duplicates (solutions within 0.1 mm)\n",
|
||
|
|
"unique_solutions = []\n",
|
||
|
|
"for sol, force in scipy_solutions:\n",
|
||
|
|
" is_duplicate = False\n",
|
||
|
|
" for unique_sol, _ in unique_solutions:\n",
|
||
|
|
" if abs(sol - unique_sol) < 0.1:\n",
|
||
|
|
" is_duplicate = True\n",
|
||
|
|
" break\n",
|
||
|
|
" if not is_duplicate:\n",
|
||
|
|
" unique_solutions.append((sol, force))\n",
|
||
|
|
"\n",
|
||
|
|
"for gap_val, force in unique_solutions:\n",
|
||
|
|
" print(f\" Gap: {gap_val:.6f} mm → Force: {force:.6f} N\")\n",
|
||
|
|
"print()"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "RLenv",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"codemirror_mode": {
|
||
|
|
"name": "ipython",
|
||
|
|
"version": 3
|
||
|
|
},
|
||
|
|
"file_extension": ".py",
|
||
|
|
"mimetype": "text/x-python",
|
||
|
|
"name": "python",
|
||
|
|
"nbconvert_exporter": "python",
|
||
|
|
"pygments_lexer": "ipython3",
|
||
|
|
"version": "3.10.9"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|