From 1cf425afe22b6c9df7176ab575dd19ad543f2ba6 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:09:02 +0100 Subject: [PATCH 01/10] relay example: update info, nicer plots --- examples/relay_feedback_system.py | 74 +++++++++++++++++++++++++------ 1 file changed, 61 insertions(+), 13 deletions(-) diff --git a/examples/relay_feedback_system.py b/examples/relay_feedback_system.py index 8de03e7..107ee46 100644 --- a/examples/relay_feedback_system.py +++ b/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() From 824c360cd2633d984f0db8ca98a3ab17a97b07aa Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:09:23 +0100 Subject: [PATCH 02/10] plot_timings: add option to export pdf --- nosnoc/plot_utils.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/nosnoc/plot_utils.py b/nosnoc/plot_utils.py index 0260e96..cb0b90d 100644 --- a/nosnoc/plot_utils.py +++ b/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() From adbfabd4e9b2856cb13378b9136bb24df4bd5f92 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:10:19 +0100 Subject: [PATCH 03/10] return empty multipliers if not there --- nosnoc/nosnoc.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nosnoc/nosnoc.py b/nosnoc/nosnoc.py index a8b4cac..cca395c 100644 --- a/nosnoc/nosnoc.py +++ b/nosnoc/nosnoc.py @@ -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] From a05bd83987eba46c3839ac5faea7e9f6d662f01a Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:10:33 +0100 Subject: [PATCH 04/10] gather multipliers in looper --- nosnoc/helpers.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/nosnoc/helpers.py b/nosnoc/helpers.py index 140ab7b..5233a63 100644 --- a/nosnoc/helpers.py +++ b/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))) From 1f9f57e7cfdea0cd889d9c6c799ebdb338dcda67 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:11:00 +0100 Subject: [PATCH 05/10] gitignore pdfs --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 28a9838..e1b8647 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,8 @@ env .vscode *.egg-info +*.pdf + private* external *~ \ No newline at end of file From 88e33484f23970896df21acf1e7517c2735e9216 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:48:36 +0100 Subject: [PATCH 06/10] add TODOs, set ipopt_tol from comp tol --- nosnoc/nosnoc.py | 10 +++++----- nosnoc/nosnoc_opts.py | 8 +++++++- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/nosnoc/nosnoc.py b/nosnoc/nosnoc.py index cca395c..f2c06b0 100644 --- a/nosnoc/nosnoc.py +++ b/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]) diff --git a/nosnoc/nosnoc_opts.py b/nosnoc/nosnoc_opts.py index c3aebd0..ea9ef80 100644 --- a/nosnoc/nosnoc_opts.py +++ b/nosnoc/nosnoc_opts.py @@ -74,7 +74,6 @@ class NosnocOpts: 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' @@ -89,7 +88,14 @@ def __repr__(self) -> str: def preprocess(self): validate(self) + self.opts_ipopt['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 + tol_ipopt = self.comp_tol * 1e-2 + self.opts_ipopt['ipopt']['tol'] = tol_ipopt + self.opts_ipopt['ipopt']['dual_inf_tol'] = tol_ipopt + self.opts_ipopt['ipopt']['compl_inf_tol'] = tol_ipopt if self.time_freezing: raise NotImplementedError() From 5b33120e9a8c47ee9e4c8a3b58e3c61c52a0bb70 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 15:49:41 +0100 Subject: [PATCH 07/10] rename: opts_ipopt -> opts_casadi_nlp --- nosnoc/nosnoc.py | 2 +- nosnoc/nosnoc_opts.py | 34 +++++++++++++++++----------------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/nosnoc/nosnoc.py b/nosnoc/nosnoc.py index f2c06b0..28bcee1 100644 --- a/nosnoc/nosnoc.py +++ b/nosnoc/nosnoc.py @@ -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=}") diff --git a/nosnoc/nosnoc_opts.py b/nosnoc/nosnoc_opts.py index ea9ef80..2792cef 100644 --- a/nosnoc/nosnoc_opts.py +++ b/nosnoc/nosnoc_opts.py @@ -63,20 +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']['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,13 +89,13 @@ 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 tol_ipopt = self.comp_tol * 1e-2 - self.opts_ipopt['ipopt']['tol'] = tol_ipopt - self.opts_ipopt['ipopt']['dual_inf_tol'] = tol_ipopt - self.opts_ipopt['ipopt']['compl_inf_tol'] = tol_ipopt + 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 if self.time_freezing: raise NotImplementedError() From 3c39e5c1adc30b9d8e4429b3faea471dcc273f8c Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 16:03:41 +0100 Subject: [PATCH 08/10] also set mu_target for IPOPT to tol_ipopt --- nosnoc/nosnoc_opts.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nosnoc/nosnoc_opts.py b/nosnoc/nosnoc_opts.py index 2792cef..7e1bdb2 100644 --- a/nosnoc/nosnoc_opts.py +++ b/nosnoc/nosnoc_opts.py @@ -96,6 +96,7 @@ def preprocess(self): 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() From f8295ea94bc3680bdc5e37831674b72589b9d225 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 16:04:37 +0100 Subject: [PATCH 09/10] add Note --- nosnoc/nosnoc_opts.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nosnoc/nosnoc_opts.py b/nosnoc/nosnoc_opts.py index 7e1bdb2..3d10a7e 100644 --- a/nosnoc/nosnoc_opts.py +++ b/nosnoc/nosnoc_opts.py @@ -92,6 +92,7 @@ def preprocess(self): 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 From 04a945a3dbe084c89af988c20ba2a727d1d5f619 Mon Sep 17 00:00:00 2001 From: Jonathan Frey Date: Fri, 9 Dec 2022 17:51:14 +0100 Subject: [PATCH 10/10] formatting, address comments --- examples/relay_feedback_system.py | 12 ++++++------ examples/simplest_example.py | 8 ++++---- nosnoc/nosnoc.py | 6 +++--- nosnoc/plot_utils.py | 1 - 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/examples/relay_feedback_system.py b/examples/relay_feedback_system.py index 107ee46..e95af5c 100644 --- a/examples/relay_feedback_system.py +++ b/examples/relay_feedback_system.py @@ -87,7 +87,6 @@ def main(): nosnoc.plot_timings(results["cpu_nlp"], figure_filename=filename) - def plot_system_3d(results): nosnoc.latexify_plot() @@ -106,6 +105,7 @@ def plot_system_3d(results): ax.grid() plt.show() + def plot_system(results): nosnoc.latexify_plot() X_sim = results["X_sim"] @@ -114,7 +114,7 @@ def plot_system(results): # state trajectory plot plt.figure() for i in range(NX): - plt.subplot(1, NX, i+1) + plt.subplot(1, NX, i + 1) plt.plot(t_grid, [x[i] for x in X_sim]) plt.grid() plt.xlabel("$t$") @@ -123,10 +123,10 @@ def plot_system(results): # 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"]))] + 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}$') diff --git a/examples/simplest_example.py b/examples/simplest_example.py index 167c32e..5fa2a29 100644 --- a/examples/simplest_example.py +++ b/examples/simplest_example.py @@ -92,10 +92,10 @@ def plot_results(results): plt.grid() # algebraic variables - 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"]))] + 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"]))] plt.subplot(3, 1, 2) n_lam = len(lambdas[0]) for i in range(n_lam): diff --git a/nosnoc/nosnoc.py b/nosnoc/nosnoc.py index 28bcee1..f8a3f9b 100644 --- a/nosnoc/nosnoc.py +++ b/nosnoc/nosnoc.py @@ -436,7 +436,7 @@ def step_equilibration(self): delta_h_ki = self.h() - self.prev_fe.h() if opts.step_equilibration == StepEquilibrationMode.HEURISTIC_MEAN: h_fe = opts.terminal_time / (opts.N_stages * opts.Nfe_list[self.ctrl_idx]) - self.cost += opts.rho_h * (self.h() -h_fe)**2 + self.cost += opts.rho_h * (self.h() - h_fe)**2 elif opts.step_equilibration == StepEquilibrationMode.HEURISTIC_DELTA: self.cost += opts.rho_h * delta_h_ki**2 elif opts.step_equilibration == StepEquilibrationMode.L2_RELAXED_SCALED: @@ -458,7 +458,7 @@ def step_equilibration(self): class NosnocSolver(NosnocFormulationObject): - # TODO: move this out of model + # TODO: move this out of solver. # NOTE: maybe dims can become part of model and are added in the preprocess function. def preprocess_model(self): # TODO: validate model @@ -666,7 +666,7 @@ def _add_primal_vector(self, symbolic: SX, lb: np.array, ub, initial): self.w0 = np.concatenate((self.w0, initial)) return - def __init__(self, opts: NosnocOpts, model: NosnocModel, ocp: Optional[NosnocOcp]=None): + def __init__(self, opts: NosnocOpts, model: NosnocModel, ocp: Optional[NosnocOcp] = None): super().__init__() diff --git a/nosnoc/plot_utils.py b/nosnoc/plot_utils.py index cb0b90d..6460059 100644 --- a/nosnoc/plot_utils.py +++ b/nosnoc/plot_utils.py @@ -45,5 +45,4 @@ def plot_timings(timings, latexify=True, title='', figure_filename=''): plt.savefig(figure_filename) print(f'stored figure as {figure_filename}') - plt.show()