Work on Relay example, add TODOs

.gitignore

@@ -3,6 +3,8 @@ env
 .vscode
 *.egg-info
 
+*.pdf
+
 private*
 external
 *~

examples/relay_feedback_system.py

@@ -1,12 +1,26 @@
 import nosnoc
 from casadi import SX, horzcat
+from datetime import datetime
+
 import numpy as np
-from mpl_toolkits import mplot3d
 import matplotlib.pyplot as plt
 
 OMEGA = 25
 XI = 0.05
 SIGMA = 1
+NX = 3
+
+## Info
+# Simulation example from
+
+# Piiroinen, Petri T., and Yuri A. Kuznetsov. "An event-driven method to simulate Filippov systems with
+# accurate computing of sliding motions." ACM Transactions on Mathematical Software (TOMS) 34.3 (2008): 1-24.
+# Equation (32)
+
+# see also:
+# M. di Bernardo, K. H. Johansson, and F. Vasca. Self-oscillations and sliding in
+# relay feedback systems: Symmetry and bifurcations. International Journal of
+# Bifurcations and Chaos, 11(4):1121-1140, 2001
 
 
 def get_relay_feedback_system_model():
@@ -45,11 +59,15 @@ def main():
     opts.cross_comp_mode = nosnoc.CrossComplementarityMode.SUM_THETAS_COMPLEMENT_WITH_EVERY_LAMBDA
     opts.step_equilibration = nosnoc.StepEquilibrationMode.HEURISTIC_MEAN
     opts.comp_tol = 1e-6
-    opts.equidistant_control_grid = False
     opts.print_level = 1
+    opts.homotopy_update_rule = nosnoc.HomotopyUpdateRule.SUPERLINEAR
+    opts.homotopy_update_exponent = 1.4
 
     Tsim = 10
     Nsim = 200
+
+    # Tsim = 1
+    # Nsim = 20
     Tstep = Tsim / Nsim
     opts.terminal_time = Tstep
 
@@ -61,11 +79,19 @@ def main():
     looper.run()
     results = looper.get_results()
 
-    plot_system(results["X_sim"], results["t_grid"])
-    nosnoc.plot_timings(results["cpu_nlp"])
+    plot_system(results)
+    # plot_system_3d(results)
+
+    filename = ""
+    filename = f"relay_timings_{datetime.utcnow().strftime('%Y-%m-%d-%H:%M:%S.%f')}.pdf"
+    nosnoc.plot_timings(results["cpu_nlp"], figure_filename=filename)
+
 
 
-def plot_system(X_sim, t_grid):
+def plot_system_3d(results):
+    nosnoc.latexify_plot()
+
+    X_sim = results["X_sim"]
     x1 = [x[0] for x in X_sim]
     x2 = [x[1] for x in X_sim]
     x3 = [x[2] for x in X_sim]
@@ -78,18 +104,40 @@ def plot_system(X_sim, t_grid):
     ax.set_ylabel("$x_2$")
     ax.set_zlabel("$x_3$")
     ax.grid()
+    plt.show()
+
+def plot_system(results):
+    nosnoc.latexify_plot()
+    X_sim = results["X_sim"]
+    t_grid = results["t_grid"]
 
     # state trajectory plot
     plt.figure()
-    plt.subplot(1, 2, 1)
-    plt.plot(t_grid, x1)
-    plt.plot(t_grid, x2)
-    plt.plot(t_grid, x3)
-    plt.xlabel("$t$")
-    plt.ylabel("$x(t)$")
+    for i in range(NX):
+        plt.subplot(1, NX, i+1)
+        plt.plot(t_grid, [x[i] for x in X_sim])
+        plt.grid()
+        plt.xlabel("$t$")
+        plt.ylabel(f"$x_{i+1}(t)$")
+
+    # algebraic variables
+    plt.figure()
+    plt.subplot(2, 1, 1)
+    lambdas = [results["lambda_sim"][0][0]
+              ] + [results["lambda_sim"][i][0] for i in range(len(results["lambda_sim"]))]
+    thetas = [results["theta_sim"][0][0]
+             ] + [results["theta_sim"][i][0] for i in range(len(results["theta_sim"]))]
+    n_lam = len(lambdas[0])
+    for i in range(n_lam):
+        plt.plot(results["t_grid"], [x[i] for x in lambdas], label=f'$\lambda_{i+1}$')
+    plt.grid()
+    plt.legend()
+
+    plt.subplot(2, 1, 2)
+    for i in range(n_lam):
+        plt.plot(results["t_grid"], [x[i] for x in thetas], label=r'$\theta_' + f'{i+1}$')
     plt.grid()
-    plt.legend(["$x_1(t)$", "$x_2(t)$", "$x_3(t)$"])
-    # TODO figure for theta/alpha
+    plt.legend()
 
     plt.show()
 

nosnoc/helpers.py

@@ -13,6 +13,7 @@ def __init__(self, solver: NosnocSolver, x0: np.ndarray, Nsim: int):
         self.time_steps = np.array([])
         self.theta_sim = []
         self.lambda_sim = []
+        self.alpha_sim = []
 
         self.cpu_nlp = np.zeros((Nsim, solver.opts.max_iter_homotopy))
 
@@ -27,6 +28,7 @@ def run(self) -> None:
             self.time_steps = np.concatenate((self.time_steps, results["time_steps"]))
             self.theta_sim += results["theta_list"]
             self.lambda_sim += results["lambda_list"]
+            self.alpha_sim += results["alpha_list"]
 
     def get_results(self) -> dict:
         self.t_grid = np.concatenate((np.array([0.0]), np.cumsum(self.time_steps)))

nosnoc/nosnoc.py

@@ -458,8 +458,10 @@ def step_equilibration(self):
 
 class NosnocSolver(NosnocFormulationObject):
 
+    # TODO: move this out of model
+    # NOTE: maybe dims can become part of model and are added in the preprocess function.
     def preprocess_model(self):
-        # Note: checks ommitted for now.
+        # TODO: validate model
         opts = self.opts
         dims = self.dims
         model = self.model
@@ -474,9 +476,7 @@ def preprocess_model(self):
         dims.n_f_sys = [model.F[i].shape[1] for i in range(dims.n_sys)]
 
         g_Stewart_list = [-model.S[i] @ model.c[i] for i in range(dims.n_sys)]
-
         g_Stewart = casadi_vertcat_list(g_Stewart_list)
-        c_all = casadi_vertcat_list(self.model.c)
 
         if opts.pss_mode == PssMode.STEP:
             # double the size of the vectors, since alpha, 1-alpha treated at same time
@@ -535,7 +535,7 @@ def preprocess_model(self):
         g_switching = SX.zeros((0, 1))
         g_convex = SX.zeros((0, 1))  # equation for the convex multiplers 1 = e' \theta
         lambda00_expr = SX.zeros(0, 0)
-        f_comp_residual = 0  # the orthogonality conditions diag(\theta) \lambda = 0.
+        f_comp_residual = SX.zeros(1)  # the orthogonality conditions diag(\theta) \lambda = 0.
 
         z = fe.rk_stage_z(0)
         if opts.pss_mode == PssMode.STEWART:
@@ -582,7 +582,7 @@ def preprocess_model(self):
         u = model.u
         model.z = z
         model.g_Stewart_fun = Function('g_Stewart_fun', [x], [g_Stewart])
-        model.c_fun = Function('c_fun', [x], [c_all])
+        model.c_fun = Function('c_fun', [x], [casadi_vertcat_list(self.model.c)])
 
         # dynamics
         model.f_x_fun = Function('f_x_fun', [x, z, u], [f_x])
@@ -755,7 +755,7 @@ def __init__(self, opts: NosnocOpts, model: NosnocModel, ocp: Optional[NosnocOcp
         # NLP Solver
         try:
             prob = {'f': self.cost, 'x': self.w, 'g': self.g, 'p': self.p}
-            self.solver = nlpsol(model.name, 'ipopt', prob, opts.opts_ipopt)
+            self.solver = nlpsol(model.name, 'ipopt', prob, opts.opts_casadi_nlp)
         except Exception as err:
             self.print_problem()
             print(f"{opts=}")
@@ -843,10 +843,10 @@ def solve(self) -> dict:
         results["lambda_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_lam]
         results["mu_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_mu]
 
-        if opts.pss_mode == PssMode.STEP:
-            results["alpha_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_alpha]
-            results["lambda_n_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_lambda_n]
-            results["lambda_p_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_lambda_p]
+        # if opts.pss_mode == PssMode.STEP:
+        results["alpha_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_alpha]
+        results["lambda_n_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_lambda_n]
+        results["lambda_p_list"] = [w_opt[flatten_layer(ind)] for ind in self.ind_lambda_p]
 
         if opts.use_fesd:
             time_steps = w_opt[self.ind_h]

nosnoc/nosnoc_opts.py

@@ -63,21 +63,20 @@ class NosnocOpts:
     equidistant_control_grid: bool = True  # NOTE: tested in test_ocp
 
     # IPOPT opts
-    opts_ipopt = dict()
-    opts_ipopt['print_time'] = 0
-    opts_ipopt['verbose'] = False
-    opts_ipopt['ipopt'] = dict()
-    opts_ipopt['ipopt']['sb'] = 'yes'
-    opts_ipopt['ipopt']['max_iter'] = 500
-    opts_ipopt['ipopt']['print_level'] = 0
-    opts_ipopt['ipopt']['bound_relax_factor'] = 0
+    opts_casadi_nlp = dict()
+    opts_casadi_nlp['print_time'] = 0
+    opts_casadi_nlp['verbose'] = False
+    opts_casadi_nlp['ipopt'] = dict()
+    opts_casadi_nlp['ipopt']['sb'] = 'yes'
+    opts_casadi_nlp['ipopt']['max_iter'] = 500
+    opts_casadi_nlp['ipopt']['print_level'] = 0
+    opts_casadi_nlp['ipopt']['bound_relax_factor'] = 0
     tol_ipopt = 1e-10
-    opts_ipopt['ipopt']['tol'] = tol_ipopt
-    opts_ipopt['ipopt']['dual_inf_tol'] = tol_ipopt
-    opts_ipopt['ipopt']['dual_inf_tol'] = tol_ipopt
-    opts_ipopt['ipopt']['compl_inf_tol'] = tol_ipopt
-    opts_ipopt['ipopt']['mu_strategy'] = 'adaptive'
-    opts_ipopt['ipopt']['mu_oracle'] = 'quality-function'
+    opts_casadi_nlp['ipopt']['tol'] = tol_ipopt
+    opts_casadi_nlp['ipopt']['dual_inf_tol'] = tol_ipopt
+    opts_casadi_nlp['ipopt']['compl_inf_tol'] = tol_ipopt
+    opts_casadi_nlp['ipopt']['mu_strategy'] = 'adaptive'
+    opts_casadi_nlp['ipopt']['mu_oracle'] = 'quality-function'
 
     time_freezing: bool = False
 
@@ -89,7 +88,16 @@ def __repr__(self) -> str:
 
     def preprocess(self):
         validate(self)
-        self.opts_ipopt['ipopt']['print_level'] = self.print_level
+
+        self.opts_casadi_nlp['ipopt']['print_level'] = self.print_level
+        # TODO: clean this up: only do this if not set by user.
+        # IPOPT tol should be smaller than outer tol, but not too much
+        # Note IPOPT option list: https://coin-or.github.io/Ipopt/OPTIONS.html
+        tol_ipopt = self.comp_tol * 1e-2
+        self.opts_casadi_nlp['ipopt']['tol'] = tol_ipopt
+        self.opts_casadi_nlp['ipopt']['dual_inf_tol'] = tol_ipopt
+        self.opts_casadi_nlp['ipopt']['compl_inf_tol'] = tol_ipopt
+        self.opts_casadi_nlp['ipopt']['mu_target'] = tol_ipopt
 
         if self.time_freezing:
             raise NotImplementedError()

nosnoc/plot_utils.py

@@ -20,7 +20,7 @@ def latexify_plot():
     matplotlib.rcParams.update(params)
 
 
-def plot_timings(timings, latexify=True, title=''):
+def plot_timings(timings, latexify=True, title='', figure_filename=''):
     # latexify plot
     if latexify:
         latexify_plot()
@@ -41,4 +41,9 @@ def plot_timings(timings, latexify=True, title=''):
     plt.xlim([-.5, Nsim - .5])
     plt.grid()
     plt.title(title)
+    if figure_filename != '':
+        plt.savefig(figure_filename)
+        print(f'stored figure as {figure_filename}')
+
+
     plt.show()