Refactor SISO uncertainty weight fit

Author: timeverettadams

example.py

@@ -0,0 +1,116 @@
+import numpy as np
+import control
+import slycot
+from matplotlib import pyplot as plt
+
+import dkpy
+
+# Frequency range
+omega_min = 0.1
+omega_max = 10
+num_omega = 100
+omega = np.logspace(np.log10(omega_min), np.log10(omega_max), num_omega)
+
+# Ground-truth system
+sys = control.tf([1, 2, 2], [1, 2.5, 1.5]) * control.tf(1, [1, 0.1])
+sys = sys * control.tf([1, 3.75, 3.5], [1, 2.5, 13])
+
+damping_ratio_1 = 0.01
+freq_nat_1 = 200
+damping_ratio_2 = 0.1
+freq_nat_2 = 500
+
+# sys = control.tf(
+#     [1 / freq_nat_1**2, 2 * damping_ratio_1 / freq_nat_1, 1],
+#     [1 / freq_nat_2**2, 2 * damping_ratio_2 / freq_nat_2, 1],
+# )
+
+# System frequency response
+frd_sys = sys.frequency_response(omega)
+fresp_sys = frd_sys.response
+mag_sys = frd_sys.magnitude
+# mag_sys = mag_sys * (1 + np.random.normal(0.0, 0.05, mag_sys.size))
+
+# Magnitude scaling
+max_mag_sys = np.max(mag_sys)
+min_mag_sys = np.min(mag_sys)
+# scale_mag_sys = np.sqrt(max_mag_sys * min_mag_sys)
+# mag_sys_nrml = mag_sys / scale_mag_sys
+
+#
+order = 5
+mosek_params = {
+    "MSK_DPAR_INTPNT_TOL_PFEAS": 1e-12,
+    "MSK_DPAR_INTPNT_TOL_DFEAS": 1e-12,
+    "MSK_DPAR_INTPNT_TOL_REL_GAP": 1e-12,
+    "MSK_DPAR_INTPNT_TOL_INFEAS": 1e-12,
+}
+linear_solver_param = {
+    "solver": "MOSEK",
+    "verbose": False,
+    "mosek_params": mosek_params,
+    "warm_start": True,
+}
+
+
+alpha = np.sqrt(omega[0] * omega[-1])
+theta = 2 * np.atan(omega / alpha)
+
+weight = dkpy.fit_magnitude_siso_ct(omega, mag_sys, 4)
+
+
+n, A, B, C, D = slycot.sb10yd(
+    discfl=0,  # Continuous-time
+    flag=1,  # Constrain stable, minimum phase
+    lendat=omega.shape[0],
+    rfrdat=np.real(fresp_sys),
+    ifrdat=np.imag(fresp_sys),
+    omega=omega,
+    n=order,
+    tol=0,  # Length of cache array
+)
+weight_slycot = control.StateSpace(A, B, C, D, dt=0)
+
+# Estimated Frequency response
+omega_fit = np.logspace(np.log10(omega_min), np.log10(omega_max), 100)
+frd_fit = weight.frequency_response(omega_fit)
+mag_fit = frd_fit.magnitude
+frd_fit_slycot = weight_slycot.frequency_response(omega_fit)
+mag_fit_slycot = frd_fit_slycot.magnitude
+
+# Plot: Preview log-magnitude response
+fig, ax = plt.subplots(layout="constrained")
+ax.semilogx(omega, control.mag2db(mag_sys), "--*", label="Magnitude Data")
+ax.semilogx(omega_fit, control.mag2db(mag_fit), label=f"Fit Magnitude (Order {order})")
+ax.semilogx(
+    omega_fit,
+    control.mag2db(mag_fit_slycot),
+    label=f"Slycot Fit Magnitude (Order {order})",
+)
+ax.set_ylabel("Magnitude (dB)")
+ax.set_xlabel("$f$ (Hz)")
+ax.legend()
+# fig.savefig(f"figs/log_cheby_fit_order_{order}.pdf")
+plt.show()
+
+# Plot: Preview magnitude response
+fig, ax = plt.subplots(layout="constrained")
+ax.semilogx(omega, mag_sys, "--*", label="Magnitude Data")
+ax.semilogx(omega_fit, mag_fit, label=f"Fit Magnitude (Order {order})")
+ax.set_ylabel("Magnitude (-)")
+ax.set_xlabel("$f$ (Hz)")
+ax.legend()
+# fig.savefig(f"figs/log_cheby_fit_order_{order}.pdf")
+plt.show()
+
+# Plot: Preview residual response
+fig, ax = plt.subplots(layout="constrained")
+ax.semilogx(omega, mag_sys**2 / mag_fit**2 - 1, "--*", label="Magnitude Data")
+ax.semilogx(
+    omega, mag_sys**2 / mag_fit_slycot**2 - 1, "--*", label=" Slycot Magnitude Data"
+)
+ax.set_ylabel("Residual (-)")
+ax.set_xlabel("$f$ (Hz)")
+ax.legend()
+# fig.savefig(f"figs/log_cheby_fit_order_{order}.pdf")
+plt.show()

src/dkpy/__init__.py

@@ -3,6 +3,7 @@
 from .controller_synthesis import *
 from .dk_iteration import *
 from .d_scale_fit import *
+from .lti_system_fit import *
 from .structured_singular_value import *
 from .uncertainty_structure import *
 from .uncertainty_characterization import *