Output switiching times (#41)

Sliding modes not handled properly yet.

examples/Acary2014/irma.py

@@ -98,13 +98,17 @@ def plot_trajectory(results, figure_filename=None):
 
     n_subplot = len(X0)
     fig, axs = plt.subplots(n_subplot, 1)
+    print(results["switch_times"])
+    xnames = ['Gal4', 'Swi5', 'Ash1', 'Cbf1', 'Gal80']
     for i in range(n_subplot):
-        axs[i].plot(results["t_grid"], results["X_sim"][:, i])
-        axs[i].hlines(thresholds[i], xmin=0, xmax=TSIM, linestyles='dotted')
+        axs[i].plot(results["t_grid"], results["X_sim"][:, i], linewidth=2)
+        axs[i].hlines(thresholds[i], xmin=0, xmax=TSIM, linestyles='dotted', linewidth=1)
         axs[i].set_xlim(0, TSIM)
         axs[i].set_ylim(0, 1.1*max(results["X_sim"][:, i]))
-        axs[i].set_ylabel(f'$x_{i+1}$')
+        axs[i].set_ylabel(xnames[i])
         axs[i].grid()
+        for t in results['switch_times']:
+            axs[i].axvline(t, linestyle="dashed", color="r", linewidth=0.5)
         if i == n_subplot - 1:
             plt.xlabel('$t$ [min]')
         else:

examples/oscilator/oscilator_example.py

@@ -72,7 +72,7 @@ def solve_oscilator(opts=None, use_g_Stewart=False, do_plot=True):
     print(f"error wrt exact solution {error:.2e}")
 
     if do_plot:
-        plot_oscilator(results["X_sim"], results["t_grid"])
+        plot_oscilator(results["X_sim"], results["t_grid"], switch_times=results["switch_times"])
     # nosnoc.plot_timings(results["cpu_nlp"])
 
     # store solution
@@ -151,16 +151,22 @@ def main_polishing():
     nosnoc.plot_timings(results["cpu_nlp"])
 
 
-def plot_oscilator(X_sim, t_grid, latexify=True):
+def plot_oscilator(X_sim, t_grid, latexify=True, switch_times=None):
     if latexify:
         nosnoc.latexify_plot()
+
+    # trajectory
     plt.figure()
     plt.subplot(1, 2, 1)
     plt.plot(t_grid, X_sim)
     plt.ylabel("$x$")
     plt.xlabel("$t$")
+    if switch_times is not None:
+        for t in switch_times:
+            plt.axvline(t, linestyle="dashed", color="k")
     plt.grid()
 
+    # state space plot
     ax = plt.subplot(1, 2, 2)
     plt.ylabel("$x_2$")
     plt.xlabel("$x_1$")

nosnoc/helpers.py

@@ -51,6 +51,7 @@ def __init__(self,
         else:
             self.print_level = solver.opts.print_level
         self.status = []
+        self.switch_times = []
 
         self.cpu_nlp = np.zeros((Nsim, solver.opts.max_iter_homotopy + 1))
 
@@ -65,6 +66,12 @@ def run(self) -> None:
                 self.solver.set("p_global", self.p_values[i, :])
             # solve
             results = self.solver.solve()
+
+            # add previous time to switch times
+            if results["switch_times"].size > 0:
+                switch_times_sim = results["switch_times"] + np.sum(self.time_steps)
+                self.switch_times += switch_times_sim.tolist()
+
             # collect
             self.X_sim += results["x_list"]
             self.xcurrent = self.X_sim[-1]
@@ -94,5 +101,6 @@ def get_results(self) -> dict:
             "w_all": self.w_all,
             "cost_vals": self.cost_vals,
             "status": self.status,
+            "switch_times": self.switch_times,
         }
         return results

nosnoc/solver.py

@@ -463,4 +463,22 @@ def get_results_from_primal_vector(prob: NosnocProblem, w_opt: np.ndarray) -> di
 
     results["v_global"] = w_opt[prob.ind_v_global]
 
+    # NOTE: this doesn't handle sliding modes well. But seems nontrivial.
+    # compute based on changes in alpha or theta
+    switch_indices = []
+    if opts.pss_mode == PssMode.STEP:
+        alpha_prev = results["alpha_list"][0]
+        for i, alpha in enumerate(results["alpha_list"][1:]):
+            if any(np.abs(alpha - alpha_prev) > 0.1):
+                switch_indices.append(i)
+            alpha_prev = alpha
+    else:
+        theta_prev = results["theta_list"][0]
+        for i, theta in enumerate(results["theta_list"][1:]):
+            if any(np.abs(theta.flatten() - theta_prev.flatten()) > 0.1):
+                switch_indices.append(i)
+            theta_prev = theta
+
+    results["switch_times"] = time_steps[switch_indices]
+
     return results