{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "594af163",
   "metadata": {},
   "source": [
    "# DCT Laboratory — Volume I, Chapter 7\n",
    "## Dynamic Enterprise Systems\n",
    "**Seed `26107`** · Companion to the chapter and AXIOM Module **AXIOM-07**\n",
    "\n",
    "Three of AXIOM-07's benches, in notebook form: the **stability spiral** of a planar\n",
    "system (spectral radius decides), the **time-semantics bench** — exact discretization\n",
    "against Euler's drift (Prop.: Exact Discretization) — and **feedback pole placement**,\n",
    "Exercise 7.12's system with the gains computed in closed form.\n",
    "Mirrored in `DCT_V1_Ch07_Lab.xlsx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f173201",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T22:20:49.343273Z",
     "iopub.status.busy": "2026-07-13T22:20:49.342989Z",
     "iopub.status.idle": "2026-07-13T22:20:49.837510Z",
     "shell.execute_reply": "2026-07-13T22:20:49.836222Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.dpi']=110\n",
    "\n",
    "import numpy as np\n",
    "SEED = 26107\n",
    "# Planar open-loop system x_{k+1} = A x_k : complex pair, stable spiral\n",
    "A = np.array([[0.60, 0.30],[-0.20, 0.90]])\n",
    "X0 = np.array([10.0, 0.0])\n",
    "\n",
    "def spectral_radius(M): return float(max(abs(np.linalg.eigvals(M))))\n",
    "def simulate(M=A, n=24, x0=X0):\n",
    "    xs=np.empty((n+1,2)); xs[0]=x0\n",
    "    for k in range(n): xs[k+1]=M@xs[k]\n",
    "    return xs\n",
    "\n",
    "# Exact vs Euler: dx/dt = lam x, sampled at DT over 5 years\n",
    "LAM, DT, N = -0.8, 0.25, 20\n",
    "def exact_path(x0=1.0):\n",
    "    return x0*np.exp(LAM*DT)**np.arange(N+1)\n",
    "def euler_path(x0=1.0):\n",
    "    return x0*(1+LAM*DT)**np.arange(N+1)\n",
    "\n",
    "# Pole placement (Ex 7.12 structure): A_c=[[0,1],[-a2,-a1]], B=[0,1]^T\n",
    "a1, a2 = 1.0, 0.5\n",
    "p1, p0 = -1.1, 0.30            # desired char poly z^2 + p1 z + p0  (roots 0.5, 0.6)\n",
    "k2 = p1 - a1                   # = -2.1\n",
    "k1 = p0 - a2                   # = -0.2\n",
    "AC = np.array([[0.0,1.0],[-a2,-a1]])\n",
    "BC = np.array([[0.0],[1.0]])\n",
    "K  = np.array([[k1,k2]])\n",
    "ACL = AC - BC@K\n",
    "\n",
    "def reference_values():\n",
    "    ex, eu = exact_path(), euler_path()\n",
    "    return {\n",
    "        \"rho_open\": round(spectral_radius(A),4),\n",
    "        \"x_t6_norm\": round(float(np.linalg.norm(simulate(n=24)[24])),4),\n",
    "        \"exact_t5\": round(float(ex[-1]),4),\n",
    "        \"euler_t5\": round(float(eu[-1]),4),\n",
    "        \"euler_drift_t5\": round(float(abs(ex[-1]-eu[-1])),4),\n",
    "        \"k1\": round(k1,4), \"k2\": round(k2,4),\n",
    "        \"rho_closed\": round(spectral_radius(ACL),4),\n",
    "    }\n",
    "if __name__ == \"__main__\":\n",
    "    [print(f\"{k:18s} {v}\") for k,v in reference_values().items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "216fdac1",
   "metadata": {},
   "source": [
    "## Panel 1 — Stability is spectral\n",
    "$x_{k+1} = A x_k$ with a complex eigenvalue pair, $\\rho(A) = 0.7746 < 1$: the\n",
    "trajectory spirals to the equilibrium at the origin. The Enterprise Stability\n",
    "Theorem in one picture — the spectrum decides, the trajectory obeys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbc05196",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T22:20:49.840223Z",
     "iopub.status.busy": "2026-07-13T22:20:49.839374Z",
     "iopub.status.idle": "2026-07-13T22:20:50.131823Z",
     "shell.execute_reply": "2026-07-13T22:20:50.130581Z"
    }
   },
   "outputs": [],
   "source": [
    "xs = simulate(n=24)\n",
    "fig, ax = plt.subplots(figsize=(6.2,5.2))\n",
    "ax.plot(xs[:,0], xs[:,1], \"o-\", c=\"#C8A24B\", lw=1.8, ms=4)\n",
    "ax.scatter([0],[0], c=\"#0B3D2E\", s=80, zorder=5, label=\"equilibrium\")\n",
    "ax.annotate(\"$x_0$\", xs[0], textcoords=\"offset points\", xytext=(8,-4))\n",
    "ax.set(xlabel=\"$x_1$\", ylabel=\"$x_2$\", title=f\"Stable spiral: ρ(A) = {spectral_radius(A):.4f} < 1\")\n",
    "ax.legend(frameon=False); ax.grid(alpha=.25); ax.set_aspect(\"equal\")\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"‖x‖ at k=24:\", round(float(np.linalg.norm(xs[24])),4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99f9be9a",
   "metadata": {},
   "source": [
    "## Panel 2 — The time-semantics bench\n",
    "$\\dot{x} = \\lambda x$ sampled at $\\Delta = 0.25$. Exact discretization multiplies by\n",
    "$e^{\\lambda\\Delta}$ each step and is **exact at the sample times**; Euler multiplies\n",
    "by $(1+\\lambda\\Delta)$ and drifts $O(\\Delta^2)$ per step. At $t=5$ the Euler path\n",
    "carries a 37% relative error — discretization is a modeling choice with a named cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f04af52c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T22:20:50.133890Z",
     "iopub.status.busy": "2026-07-13T22:20:50.133700Z",
     "iopub.status.idle": "2026-07-13T22:20:50.336004Z",
     "shell.execute_reply": "2026-07-13T22:20:50.334513Z"
    }
   },
   "outputs": [],
   "source": [
    "t = np.arange(N+1)*DT\n",
    "ex, eu = exact_path(), euler_path()\n",
    "fig, ax = plt.subplots(figsize=(8,4.2))\n",
    "ax.plot(t, ex, c=\"#0B3D2E\", lw=2.4, label=\"exact discretization  $e^{\\\\lambda\\\\Delta}$\")\n",
    "ax.plot(t, eu, c=\"#C8A24B\", lw=2.2, ls=\"--\", label=\"Euler  $(1+\\\\lambda\\\\Delta)$\")\n",
    "ax.set(xlabel=\"years\", ylabel=\"x\", title=\"Same continuous law, two discrete semantics (seed 26107)\")\n",
    "ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"exact t=5: {ex[-1]:.4f}   euler t=5: {eu[-1]:.4f}   drift: {abs(ex[-1]-eu[-1]):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f385271c",
   "metadata": {},
   "source": [
    "## Panel 3 — Feedback places the poles\n",
    "Exercise 7.12's system: $A_c = \\begin{bmatrix}0&1\\\\-a_2&-a_1\\end{bmatrix}$,\n",
    "$B = [0,1]^\\top$. Desired eigenvalues $\\{0.5, 0.6\\}$; in controllable canonical\n",
    "form the gains read off the characteristic coefficients: $k_1 = p_0 - a_2$,\n",
    "$k_2 = p_1 - a_1$. The Feedback Stability Theorem, constructively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af6510fd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T22:20:50.338166Z",
     "iopub.status.busy": "2026-07-13T22:20:50.337931Z",
     "iopub.status.idle": "2026-07-13T22:20:50.575034Z",
     "shell.execute_reply": "2026-07-13T22:20:50.573337Z"
    }
   },
   "outputs": [],
   "source": [
    "print(\"K =\", K.ravel(), \"  closed-loop eigenvalues:\", np.round(np.linalg.eigvals(ACL),4))\n",
    "print(\"ρ(A−BK) =\", round(spectral_radius(ACL),4))\n",
    "xs_o = simulate(AC, 20, np.array([1.0,0.0]))\n",
    "xs_c = simulate(ACL, 20, np.array([1.0,0.0]))\n",
    "fig, ax = plt.subplots(figsize=(8,4.0))\n",
    "ax.plot(np.linalg.norm(xs_o,axis=1), c=\"#8A8F8B\", lw=2, label=\"open loop ‖x‖\")\n",
    "ax.plot(np.linalg.norm(xs_c,axis=1), c=\"#C8A24B\", lw=2.3, label=\"closed loop ‖x‖ (poles at 0.5, 0.6)\")\n",
    "ax.set(xlabel=\"k\", ylabel=\"‖x‖\", title=\"Feedback modifies trajectories — and their decay rate\")\n",
    "ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad540bd6",
   "metadata": {},
   "source": [
    "## Validation — agrees with `DCT_V1_Ch07_Lab.xlsx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a529be3f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T22:20:50.577107Z",
     "iopub.status.busy": "2026-07-13T22:20:50.576898Z",
     "iopub.status.idle": "2026-07-13T22:20:50.584286Z",
     "shell.execute_reply": "2026-07-13T22:20:50.583215Z"
    }
   },
   "outputs": [],
   "source": [
    "ref = reference_values()\n",
    "expected = {\"rho_open\":0.7746,\"x_t6_norm\":0.0254,\"exact_t5\":0.0183,\"euler_t5\":0.0115,\n",
    " \"euler_drift_t5\":0.0068,\"k1\":-0.2,\"k2\":-2.1,\"rho_closed\":0.6}\n",
    "for k,v in expected.items():\n",
    "    assert abs(ref[k]-v)<5e-4, f\"MISMATCH {k}\"\n",
    "    print(f\"PASS  {k:18s} {ref[k]}\")\n",
    "print(\"\\nAll checkpoints agree — seed 26107.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce118b6c",
   "metadata": {},
   "source": [
    "**Next**: Exercises 7.9–7.12 (★ exercises live here — the rotation matrix and the pole placement); AXIOM-07's cobweb theater and time-semantics bench extend the panels. Solutions: IM Ch. 7."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}