Skip to content

Commit

Permalink
Order plots and minor fixes (#42)
Browse files Browse the repository at this point in the history
  • Loading branch information
@FreyJo
FreyJo authored Mar 23, 2023
1 parent 76c3258 commit 6bf9206
Show file tree
Hide file tree
Showing 6 changed files with 221 additions and 50 deletions.
78 changes: 49 additions & 29 deletions examples/Acary2014/irma_integration_order_experiment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
NFE_VALUES = [3]
# NSIM_VALUES = [1, 3, 10, 20, 50, 100, 300] # convergence issues for Legendre
NSIM_VALUES = [1, 3, 9, 18, 50, 100, 300]
USE_FESD_VALUES = [True, False]
# USE_FESD_VALUES = [False]


# # NOTE: this has convergence issues
Expand All @@ -27,6 +29,10 @@


def pickle_results(results, filename):
# create directory if it does not exist
if not os.path.exists(BENCHMARK_DATA_PATH):
os.makedirs(BENCHMARK_DATA_PATH)
# save
file = os.path.join(BENCHMARK_DATA_PATH, filename)
with open(file, 'wb') as f:
pickle.dump(results, f)
Expand Down Expand Up @@ -71,18 +77,21 @@ def get_results_filename(opts):
filename += 'dt' + str(opts.terminal_time) + '_'
filename += 'Tsim' + str(TSIM) + '_'
filename += opts.irk_scheme.name
if not opts.use_fesd:
filename += '_nofesd'
filename += '.pickle'
return filename


def get_opts(Nsim, n_s, N_fe, scheme):
def get_opts(Nsim, n_s, N_fe, scheme, use_fesd):
opts = get_default_options()
Tstep = TSIM / Nsim
opts.terminal_time = Tstep
opts.comp_tol = TOL
opts.sigma_N = TOL
opts.irk_scheme = scheme
opts.print_level = 1
opts.use_fesd = use_fesd

opts.n_s = n_s
opts.N_finite_elements = N_fe
Expand All @@ -92,10 +101,10 @@ def get_opts(Nsim, n_s, N_fe, scheme):


def run_benchmark():
for n_s, N_fe, Nsim, scheme in itertools.product(NS_VALUES, NFE_VALUES, NSIM_VALUES, SCHEMES):
for n_s, N_fe, Nsim, scheme, use_fesd in itertools.product(NS_VALUES, NFE_VALUES, NSIM_VALUES, SCHEMES, USE_FESD_VALUES):

model = get_irma_model(SWITCH_ON, LIFTING)
opts = get_opts(Nsim, n_s, N_fe, scheme)
opts = get_opts(Nsim, n_s, N_fe, scheme, use_fesd)
solver = nosnoc.NosnocSolver(opts, model)

# loop
Expand Down Expand Up @@ -128,41 +137,52 @@ def order_plot():
ref_results = get_reference_solution()
x_end_ref = ref_results['X_sim'][-1]

linestyles = ['-', '--', '-.', ':']
linestyles = ['-', '--', '-.', ':', ':', '-.', '--', '-']
marker_types = ['o', 's', 'v', '^', '>', '<', 'd', 'p']
# SCHEME = nosnoc.IrkSchemes.RADAU_IIA
SCHEME = nosnoc.IrkSchemes.RADAU_IIA

plt.figure()
for i, n_s in enumerate(NS_VALUES):
errors = []
step_sizes = []
for Nsim in NSIM_VALUES:
opts = get_opts(Nsim, n_s, N_fe, SCHEME)
filename = get_results_filename(opts)
results = unpickle_results(filename)
x_end = results['X_sim'][-1]
n_fail = count_failures(results)
error = np.max(np.abs(x_end - x_end_ref))
print("opts.n_s: ", opts.n_s, "opts.terminal_time: ", opts.terminal_time, "error: ", error, "n_fail: ", n_fail)
errors.append(error)
step_sizes.append(opts.terminal_time)

label = r'$n_s=' + str(n_s) +'$'
if results['opts'].irk_scheme == nosnoc.IrkSchemes.RADAU_IIA:
if n_s == 1:
label = 'implicit Euler: 1'
ax = plt.figure()
for use_fesd in [True, False]:
for i, n_s in enumerate(NS_VALUES):
errors = []
step_sizes = []
for Nsim in NSIM_VALUES:
opts = get_opts(Nsim, n_s, N_fe, SCHEME, use_fesd)
filename = get_results_filename(opts)
results = unpickle_results(filename)
print(f"loading filde {filename}")
x_end = results['X_sim'][-1]
n_fail = count_failures(results)
error = np.max(np.abs(x_end - x_end_ref))
print("opts.n_s: ", opts.n_s, "opts.terminal_time: ", opts.terminal_time, "error: ", error, "n_fail: ", n_fail)
errors.append(error)
step_sizes.append(opts.terminal_time)

label = r'$n_s=' + str(n_s) +'$'
if results['opts'].irk_scheme == nosnoc.IrkSchemes.RADAU_IIA:
if n_s == 1:
label = 'implicit Euler: 1'
else:
label = 'Radau IIA: ' + str(2*n_s-1)
elif results['opts'].irk_scheme == nosnoc.IrkSchemes.GAUSS_LEGENDRE:
label = 'Gauss-Legendre: ' + str(2*n_s)
if use_fesd:
label += ', FESD'
else:
label = 'Radau IIA: ' + str(2*n_s-1)
elif results['opts'].irk_scheme == nosnoc.IrkSchemes.GAUSS_LEGENDRE:
label = 'Gauss-Legendre: ' + str(2*n_s)
plt.plot(step_sizes, errors, label=label, marker='o', linestyle=linestyles[i])
label += ', Standard'
plt.plot(step_sizes, errors, label=label, marker=marker_types[i], linestyle=linestyles[i])
plt.grid()
plt.xlabel('Step size')
plt.ylabel('Error')
plt.yscale('log')
plt.xscale('log')
plt.legend()
# plt.legend(loc='center left')
ax.legend(loc='upper left', bbox_to_anchor=(0.05, .97), ncol=2, framealpha=1.0)

plt.savefig(f'irma_benchmark_{SCHEME.name}.pdf', bbox_inches='tight')
fig_filename = f'irma_benchmark_{SCHEME.name}.pdf'
plt.savefig(fig_filename, bbox_inches='tight')
print(f"Saved figure to {fig_filename}")
plt.show()


Expand Down
151 changes: 151 additions & 0 deletions examples/oscillator/integration_order_experiment_oscillator.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
from matplotlib import pyplot as plt
import numpy as np
from oscillator_example import get_oscillator_model, X_SOL, TSIM
import nosnoc
import pickle
import os
import itertools

BENCHMARK_DATA_PATH = 'private_oscillator_benchmark_data'
SCHEMES = [nosnoc.IrkSchemes.GAUSS_LEGENDRE, nosnoc.IrkSchemes.RADAU_IIA]

NS_VALUES = [1, 2, 3, 4]
NFE_VALUES = [2]
NSIM_VALUES = [1, 5, 9, 10, 20, 40, 60, 100]
USE_FESD_VALUES = [True, False]

TOL = 1e-12


def pickle_results(results, filename):
# create directory if it does not exist
if not os.path.exists(BENCHMARK_DATA_PATH):
os.makedirs(BENCHMARK_DATA_PATH)
# save
file = os.path.join(BENCHMARK_DATA_PATH, filename)
with open(file, 'wb') as f:
pickle.dump(results, f)

def unpickle_results(filename):
file = os.path.join(BENCHMARK_DATA_PATH, filename)
with open(file, 'rb') as f:
results = pickle.load(f)
return results


def get_results_filename(opts: nosnoc.NosnocOpts):
filename = 'oscillator_bm_results_'
filename += 'Nfe_' + str(opts.N_finite_elements) + '_'
filename += 'ns' + str(opts.n_s) + '_'
filename += 'tol' + str(opts.comp_tol) + '_'
filename += 'dt' + str(opts.terminal_time) + '_'
filename += 'Tsim' + str(TSIM) + '_'
filename += opts.irk_scheme.name + '_'
filename += opts.pss_mode.name
if not opts.use_fesd:
filename += '_nofesd'
filename += '.pickle'
return filename


def get_opts(Nsim, n_s, N_fe, scheme, use_fesd):
opts = nosnoc.NosnocOpts()
Tstep = TSIM / Nsim
opts.terminal_time = Tstep
opts.comp_tol = TOL
opts.sigma_N = TOL
opts.irk_scheme = scheme
opts.print_level = 0
opts.tol_ipopt = TOL

opts.n_s = n_s
opts.N_finite_elements = N_fe
opts.pss_mode = nosnoc.PssMode.STEP
opts.use_fesd = use_fesd
return opts


def run_benchmark():
for n_s, N_fe, Nsim, scheme, use_fesd in itertools.product(NS_VALUES, NFE_VALUES, NSIM_VALUES, SCHEMES, USE_FESD_VALUES):

model = get_oscillator_model()
opts = get_opts(Nsim, n_s, N_fe, scheme, use_fesd)
solver = nosnoc.NosnocSolver(opts, model)

# loop
looper = nosnoc.NosnocSimLooper(solver, model.x0, Nsim, print_level=1)
looper.run()
results = looper.get_results()
results['opts'] = opts
filename = get_results_filename(opts)
pickle_results(results, filename)
del solver, looper, results, model, opts



def count_failures(results):
status_list: list = results['status']
return len([x for x in status_list if x != nosnoc.Status.SUCCESS])



def order_plot():
nosnoc.latexify_plot()
N_fe = 2

linestyles = ['-', '--', '-.', ':', ':', '-.', '--', '-']
marker_types = ['o', 's', 'v', '^', '>', '<', 'd', 'p']
SCHEME = nosnoc.IrkSchemes.GAUSS_LEGENDRE
# SCHEME = nosnoc.IrkSchemes.RADAU_IIA

ax = plt.figure()
for use_fesd in [True, False]:
for i, n_s in enumerate(NS_VALUES):
errors = []
step_sizes = []
for Nsim in NSIM_VALUES:
opts = get_opts(Nsim, n_s, N_fe, SCHEME, use_fesd)
filename = get_results_filename(opts)
results = unpickle_results(filename)
print(f"loading filde {filename}")
x_end = results['X_sim'][-1]
n_fail = count_failures(results)
error = np.max(np.abs(x_end - X_SOL))
print("opts.n_s: ", opts.n_s, "opts.terminal_time: ", opts.terminal_time, "error: ", error, "n_fail: ", n_fail)
errors.append(error)
step_sizes.append(opts.terminal_time)

label = r'$n_s=' + str(n_s) +'$'
if results['opts'].irk_scheme == nosnoc.IrkSchemes.RADAU_IIA:
if n_s == 1:
label = 'implicit Euler: 1'
else:
label = 'Radau IIA: ' + str(2*n_s-1)
elif results['opts'].irk_scheme == nosnoc.IrkSchemes.GAUSS_LEGENDRE:
label = 'Gauss-Legendre: ' + str(2*n_s)
if use_fesd:
label += ', FESD'
else:
label += ', Standard'
plt.plot(step_sizes, errors, label=label, marker=marker_types[i], linestyle=linestyles[i])
plt.grid()
plt.xlabel('Step size')
plt.ylabel('Error')
plt.yscale('log')
plt.xscale('log')
# plt.legend(loc='center left')
plt.legend()
# ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.0), ncol=2, framealpha=1.0)

fig_filename = f'oscillator_benchmark_{SCHEME.name}.pdf'
plt.savefig(fig_filename, bbox_inches='tight')
print(f"Saved figure to {fig_filename}")
plt.show()


if __name__ == "__main__":
# generate data
# run_benchmark()

# evalute
order_plot()
22 changes: 11 additions & 11 deletions examples/oscilator/oscilator_example.py → examples/oscillator/oscillator_example.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
TSIM = np.pi / 2
X_SOL = np.array([np.exp(TSIM-1) * np.cos(2*np.pi * (TSIM-1)), -np.exp(TSIM-1) * np.sin(2*np.pi*(TSIM-1))])

def get_oscilator_model(use_g_Stewart=False):
def get_oscillator_model(use_g_Stewart=False):

# Initial Value
x0 = np.array([np.exp([-1])[0], 0])
Expand Down Expand Up @@ -51,11 +51,11 @@ def get_default_options():
return opts


def solve_oscilator(opts=None, use_g_Stewart=False, do_plot=True):
def solve_oscillator(opts=None, use_g_Stewart=False, do_plot=True):
if opts is None:
opts = get_default_options()

model = get_oscilator_model(use_g_Stewart)
model = get_oscillator_model(use_g_Stewart)

Nsim = 29
Tstep = TSIM / Nsim
Expand All @@ -72,12 +72,12 @@ 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"], switch_times=results["switch_times"])
plot_oscillator(results["X_sim"], results["t_grid"], switch_times=results["switch_times"])
# nosnoc.plot_timings(results["cpu_nlp"])

# store solution
# import json
# json_file = 'oscilator_results_ref.json'
# json_file = 'oscillator_results_ref.json'
# with open(json_file, 'w') as f:
# json.dump(results['w_sim'], f, indent=4, sort_keys=True, default=make_object_json_dumpable)
# print(f"saved results in {json_file}")
Expand All @@ -87,7 +87,7 @@ def main_least_squares():

# load reference solution
# import json
# json_file = 'oscilator_results_ref.json'
# json_file = 'oscillator_results_ref.json'
# with open(json_file, 'r') as f:
# w_sim_ref = json.load(f)

Expand All @@ -109,7 +109,7 @@ def main_least_squares():
# opts.homotopy_update_rule = nosnoc.HomotopyUpdateRule.SUPERLINEAR
opts.homotopy_update_slope = 0.1

model = get_oscilator_model()
model = get_oscillator_model()

Nsim = 29
Tstep = TSIM / Nsim
Expand All @@ -129,7 +129,7 @@ def main_least_squares():
print(f"error wrt exact solution {error:.2e}")

breakpoint()
plot_oscilator(results["X_sim"], results["t_grid"])
plot_oscillator(results["X_sim"], results["t_grid"])
# nosnoc.plot_timings(results["cpu_nlp"])


Expand All @@ -145,13 +145,13 @@ def main_polishing():
opts.mpcc_mode = nosnoc.MpccMode.FISCHER_BURMEISTER
opts.print_level = 3

results = solve_oscilator(opts, do_plots=False)
results = solve_oscillator(opts, do_plots=False)
print(f"max cost_val = {max(results['cost_vals']):.2e}")

nosnoc.plot_timings(results["cpu_nlp"])


def plot_oscilator(X_sim, t_grid, latexify=True, switch_times=None):
def plot_oscillator(X_sim, t_grid, latexify=True, switch_times=None):
if latexify:
nosnoc.latexify_plot()

Expand Down Expand Up @@ -198,6 +198,6 @@ def make_object_json_dumpable(input):
raise TypeError(f"Cannot make input of type {type(input)} dumpable.")

if __name__ == "__main__":
solve_oscilator(use_g_Stewart=False, do_plot=True)
solve_oscillator(use_g_Stewart=False, do_plot=True)
# main_least_squares()
# main_polishing()
2 changes: 1 addition & 1 deletion nosnoc/problem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -483,7 +483,7 @@ def create_complementarity_constraints(self, sigma_p: ca.SX, tau: ca.SX, Uk: ca.
if not opts.use_fesd:
for j in range(opts.n_s):
self.create_complementarity([self.Lambda(stage=j)], self.Theta(stage=j), sigma_p,
tau, Uk)
tau)
elif opts.cross_comp_mode == CrossComplementarityMode.COMPLEMENT_ALL_STAGE_VALUES_WITH_EACH_OTHER:
for j in range(opts.n_s):
# cross comp with prev_fe
Expand Down
2 changes: 1 addition & 1 deletion nosnoc/solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -479,6 +479,6 @@ def get_results_from_primal_vector(prob: NosnocProblem, w_opt: np.ndarray) -> di
switch_indices.append(i)
theta_prev = theta

results["switch_times"] = time_steps[switch_indices]
results["switch_times"] = np.array([time_steps[i] for i in switch_indices])

return results
Loading

0 comments on commit 6bf9206

Please sign in to comment.