From 5acf66ab558a9b97a0ff427e141eb681ccced129 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Wed, 25 Sep 2024 07:53:18 +0200
Subject: [PATCH 01/10] variable measurements and repeats

---
 ...ontiguous_parameter_identification_wang.py | 323 ++++++++++++++++++
 1 file changed, 323 insertions(+)
 create mode 100644 examples-gallery/plot_non_contiguous_parameter_identification_wang.py

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
new file mode 100644
index 00000000..f3341232
--- /dev/null
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -0,0 +1,323 @@
+# %%
+"""
+Parameter Identification from Noncontiguous Measurements with Bias
+==================================================================
+
+In parameter estimation it is common to collect measurements of a system's
+trajectories from distinct experiments. For example, if you are identifying the
+parameters of a mass-spring-damper system, you may excite the system with
+different initial conditions multiple times and measure the position of the
+mass. The data cannot simply be stacked end-to-end in time and the
+identification run because the measurement data would be discontinuous between
+each measurement trial.
+
+A workaround in opty is to create a set of differential equations with unique
+state variables for each measurement trial that all share the same constant
+parameters. You can then identify the parameters from all measurement trials
+simultaneously by passing the uncoupled differential equations to opty.
+
+Mass-spring-damper-friction Example
+===================================
+
+The position of a simple system consisting of a mass connected to a fixed point
+by a spring and a damper and subject to friction is simulated and recorded as
+noisy measurements with bias. The spring constant, the damping coefficient
+and the coefficient of friction will be identified.
+
+**State Variables**
+
+- :math:`x_i`: position of the mass of the i-th measurement trial [m]
+- :math:`u_i`: speed of the mass of the i-th measurement trial [m/s]
+
+**Parameters**
+
+- :math:`m`: mass [kg]
+- :math:`c`: linear damping coefficient [Ns/m]
+- :math:`k`: linear spring constant [N/m]
+- :math:`l_0`: natural length of the spring [m]
+- :math:`friction`: friction coefficient [N]
+
+"""
+import sympy as sm
+import sympy.physics.mechanics as me
+import numpy as np
+from scipy.integrate import solve_ivp
+from opty import Problem
+from opty.utils import parse_free
+import matplotlib.pyplot as plt
+
+# %%
+# Equations of Motion.
+# --------------------
+#
+number_of_measurements = 4
+number_of_repeats = 3
+t0, tf = 0.0, 10.0
+num_nodes = 500
+
+m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction')
+bias_ident = sm.symbols(f'bias:{number_of_measurements}')
+
+t = me.dynamicsymbols._t
+x = sm.Matrix([me.dynamicsymbols([f'x{i}_{j}' for j in range(number_of_repeats)])
+        for i in range(number_of_measurements)])
+u = sm.Matrix([me.dynamicsymbols([f'u{i}_{j}' for j in range(number_of_repeats)])
+        for i in range(number_of_measurements)])
+
+n = sm.Matrix([me.dynamicsymbols([f'n{i}_{j}' for j in range(number_of_repeats)])
+        for i in range(number_of_measurements)])
+
+xh, uh = me.dynamicsymbols('xh, uh')
+
+# %%
+# Form the equations of motion, such that the first number_of_repeats equations
+# belong to the first measurement, the second number_of_repeats equations belong
+# to the second measurement, and so on.
+def kd_eom_rh(xh, uh, m, c, k, l0, friction):
+    """" sets up the eoms for the system"""
+    N = me.ReferenceFrame('N')
+    O, P = sm.symbols('O, P', cls=me.Point)
+    O.set_vel(N, 0)
+
+    P.set_pos(O, xh*N.x)
+    P.set_vel(N, uh*N.x)
+    bodies = [me.Particle('part', P, m)]
+    # :math:`\tanh(\alpha \cdot x) \approx sign(x)` for large :math:`\alpha`, and
+    # is differentiale everywhere.
+    forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(30*uh)*N.x)]
+    kd = sm.Matrix([uh - xh.diff(t)])
+    KM = me.KanesMethod(N,
+                        q_ind=[xh],
+                        u_ind=[uh],
+                        kd_eqs=kd
+    )
+    fr, frstar = KM.kanes_equations(bodies, forces)
+    rhs = KM.rhs()
+    return (kd, fr + frstar, rhs)
+
+kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction)
+
+kd_total = sm.Matrix([])
+eom_total = sm.Matrix([])
+rhs_total = sm.Matrix([])
+for i in range(number_of_measurements):
+    for j in range(number_of_repeats):
+        eom_h = me.msubs(fr_frstar, {xh: x[i, j], uh: u[i, j], uh.diff(t): u[i, j].diff(t)})
+        eom_h = eom_h + sm.Matrix([n[i, j]])
+        kd_h = kd.subs({xh: x[i, j], uh: u[i, j]})
+        kd_total = kd_total.col_join(kd_h)
+        eom_total = eom_total.col_join(eom_h)
+        rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})])
+        rhs_total = rhs_total.col_join(rhs_h)
+
+eom = kd_total.col_join(eom_total)
+uh =sm.Matrix([u[i, j] for i in range(number_of_measurements)
+    for j in range(number_of_repeats)])
+rhs = uh.col_join(rhs_total)
+
+for i in range(number_of_repeats):
+    sm.pprint(eom[i])
+for i in range(3):
+    print('.........')
+for i in range(number_of_repeats):
+    sm.pprint(eom[-(number_of_repeats-i)])
+
+# %%
+# Generate Noisy Measurement Data with Bias
+# -----------------------------------------
+# Create 'number_of_measurements' sets of measurements with different initial
+# conditions. To get the measurements, the equations of motion are integrated,
+# and then noise is added to each point in time of the solution. Finally a bias
+# is added.
+#
+
+states = [x[i, j] for i in range(number_of_measurements)
+                for j in range(number_of_repeats)] + \
+        [u[i, j] for i in range(number_of_measurements)
+                for j in range(number_of_repeats)]
+parameters = [m, c, k, l0, friction]
+par_vals = [1.0, 0.25, 2.0, 1.0, 0.5]
+
+eval_rhs = sm.lambdify(states + parameters, rhs)
+
+times = np.linspace(t0, tf, num=num_nodes)
+
+measurements = []
+# %%
+
+for i in range(number_of_measurements):
+    for j in range(number_of_repeats):
+        seed = 1234*(i+1)*(j+1)
+        np.random.seed(seed)
+        x0 = 4*np.random.randn(2*number_of_measurements * number_of_repeats)
+        sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(),
+                    (t0, tf), x0, t_eval=times, args=(par_vals,))
+    seed = 10000 + 12345*i
+    np.random.seed(seed)
+    measurements.append(sol.y[1, :] + 1.0*np.random.randn(num_nodes))
+measurements = np.array(measurements)
+print(measurements.shape)
+print(sol.y.shape)
+
+fig, ax = plt.subplots(2*number_of_measurements*number_of_repeats, 1,
+    figsize=(6.4, 1.5*number_of_measurements*number_of_repeats),
+    sharex=True, constrained_layout=True)
+
+for i in range(2*number_of_measurements*number_of_repeats):
+    ax[i].plot(times, sol.y[i, :], color='blue')
+    ax[i].set_ylabel(f'{states[i]}')
+
+# Add bias to the measurements
+bias = 0.0 *np.random.uniform(0., 1., size=number_of_measurements)
+measurements = np.array(measurements) + bias[:, np.newaxis]
+
+print(measurements.shape)
+
+# %%
+# Setup the Identification Problem
+# --------------------------------
+#
+# The goal is to identify the damping coefficient :math:`c`, the spring
+# constant :math:`k` and the coefficient of friction :math:`friction`.
+# The objective :math:`J` is to minimize the least square
+# difference in the optimal simulation as compared to the measurements.  If
+# some measurement is considered more reliable, its weight :math:`w` may be
+# increased relative to the other measurements.
+#
+#
+# objective = :math:`\int_{t_0}^{t_f} \left[ \sum_{i=1}^{i = \text{number_of_measurements}} (w_i (x_i - x_i^m)^2 \right] \hspace{2pt} dt`
+#
+interval_value = (tf - t0) / (num_nodes - 1)
+
+w =[1 for _ in range(number_of_measurements)]
+
+nm = number_of_measurements
+nr = number_of_repeats
+def obj(free):
+    sum = 0.0
+    for i in range(nm):
+        for j in range(nr):
+            sum += (w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] -
+                measurements[i])**2)
+
+
+def obj_grad(free):
+    grad = np.zeros_like(free)
+    for i in range(nm):
+        for j in range(nr):
+            grad[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] = 2*w[i]*interval_value*(
+                free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - measurements[i]
+            )
+    return grad
+
+
+# %%
+# By not including :math:`c`, :math:`k`, :math:`friction`, :math:`bias_i`
+# in the parameter map, they will be treated as unknown parameters.
+par_map = {m: par_vals[0], l0: par_vals[3]}
+
+# %%
+bounds = {
+    c: (0.01, 2.0),
+    k: (0.1, 10.0),
+    friction: (0.1, 10.0),
+}
+noise_scale = 2.0
+known_trajectory_map = {}
+for i in range(number_of_measurements):
+    for j in range(number_of_repeats):
+        seed = 10000 + 12345*(i+1)*(j+1)
+        np.random.seed(seed)
+        known_trajectory_map = (known_trajectory_map |
+            {n[i, j]: noise_scale*np.random.randn(num_nodes)})
+print(len(known_trajectory_map))
+print('states', len(states))
+problem = Problem(
+    obj,
+    obj_grad,
+    eom,
+    states,
+    num_nodes,
+    interval_value,
+    known_parameter_map=par_map,
+    time_symbol=me.dynamicsymbols._t,
+    integration_method='midpoint',
+    bounds=bounds,
+    known_trajectory_map=known_trajectory_map
+)
+print(problem.num_free)
+# %%
+# This give the sequence of the unknown parameters at the tail of solution.
+print(problem.collocator.unknown_parameters)
+print(problem.collocator.unknown_input_trajectories)
+# %%
+# Create an Initial Guess
+# -----------------------
+#
+# It is reasonable to use the measurements as initial guess for the states
+# because they would be available. Here, only the measurements of the position
+# are used and the speeds are set to zero. The last values are the guesses
+# for :math:`bias_i`, :math:`c` and :math:`k`, respectively.
+#
+initial_guess = np.array([])
+for i in range(number_of_measurements):
+    for j in range(number_of_repeats):
+        initial_guess = np.concatenate((initial_guess, measurements[i, :]))
+initial_guess = np.hstack((initial_guess,
+                np.zeros(number_of_measurements*number_of_repeats*num_nodes),
+                [0.1, 3.0, 2.0]))
+
+print(len(initial_guess))
+
+# %%
+# Solve the Optimization Problem
+# ------------------------------
+#
+# Here, the initial guess is the result of some previous run, stored in
+# 'non_contiguous_parameter_identification_bias_solution.npy'.
+#initial_guess = np.load('non_contiguous_parameter_identification_bias_solution.npy')
+solution, info = problem.solve(initial_guess)
+# %%
+# This is how the solution may be saved for a future initial guess
+# ```np.save('non_contiguous_parameter_identification_bias_solution', solution)```
+
+print(info['status_msg'])
+print(f'final value of the objective function is {info['obj_val']:.2f}' )
+
+# %%
+problem.plot_objective_value()
+
+# %%
+problem.plot_constraint_violations(solution)
+
+
+# %%
+# The identified parameters are:
+# ------------------------------
+#
+print(f'Estimate of damping coefficient is      {solution[-3]: 1.2f}' +
+      f' Percentage error is {(solution[-3]-par_vals[1])/solution[-3]*100:1.2f} %')
+print(f'Estimate of the spring constant is      {solution[-1]: 1.2f}' +
+      f' Percentage error is {(solution[-1]-par_vals[2])/solution[-1]*100:1.2f} %')
+print(f'Estimate of the friction coefficient is {solution[-2]: 1.2f}'
+      f' Percentage error is {(solution[-2]-par_vals[-1])/solution[-2]*100:1.2f} %')
+
+# %%
+# Plot the measurements and the trajectories calculated.
+#
+fig, ax = plt. subplots(number_of_measurements, 1,
+    figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True)
+for i in range(number_of_measurements):
+    ax[i].plot(times, solution[i*num_nodes:(i+1)*num_nodes], color='red')
+    ax[i].plot(times, measurements[i], color='blue', lw=0.5)
+    ax[i].set_ylabel(f'{states[i]}')
+ax[-1].set_xlabel('Time [s]')
+ax[0].set_title('Trajectories')
+prevent_output = True
+
+# %%
+#
+# sphinx_gallery_thumbnail_number = 3
+problem.plot_trajectories(solution)
+
+plt.show()

From 231908bcb6cd7555168e2030f73c84c31494b572 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Wed, 25 Sep 2024 13:15:11 +0200
Subject: [PATCH 02/10] including noise for every measurement

---
 ...ontiguous_parameter_identification_wang.py | 136 ++++++++++--------
 1 file changed, 75 insertions(+), 61 deletions(-)

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
index f3341232..8aeb22b7 100644
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -1,7 +1,6 @@
-# %%
 """
-Parameter Identification from Noncontiguous Measurements with Bias
-==================================================================
+Parameter Identification from Noncontiguous Measurements with Bias and added Noise
+==================================================================================
 
 In parameter estimation it is common to collect measurements of a system's
 trajectories from distinct experiments. For example, if you are identifying the
@@ -15,6 +14,18 @@
 state variables for each measurement trial that all share the same constant
 parameters. You can then identify the parameters from all measurement trials
 simultaneously by passing the uncoupled differential equations to opty.
+This idea is due to Jason Moore.
+
+In addition, for each measurement, one may create a set of differential
+equations with unique state variables that again share the same constant
+parameters. In addition, one adds noise to these equations. (Intuitively this
+noise corresponds to a force acting on the system).
+This idea is due to Huawei Wang and A. J. van den Bogert.
+
+So, if one has 'number_of_measurements' measurements, and one creates
+'number_of_repeats' differential equations for each measurement, one will have
+'number_of_measurements' * 'number_of_repeats' uncoupled differential equations,
+all sharing the same constant parameters.
 
 Mass-spring-damper-friction Example
 ===================================
@@ -26,8 +37,12 @@
 
 **State Variables**
 
-- :math:`x_i`: position of the mass of the i-th measurement trial [m]
-- :math:`u_i`: speed of the mass of the i-th measurement trial [m/s]
+- :math:`x_{i, j}`: position of the mass of the i-th measurement trial [m]
+- :math:`u_{i, j}`: speed of the mass of the i-th measurement trial [m/s]
+
+**Noise Variables**
+
+- :math:`n_{i, j}`: noise added [N]
 
 **Parameters**
 
@@ -49,14 +64,15 @@
 # %%
 # Equations of Motion.
 # --------------------
-#
-number_of_measurements = 4
-number_of_repeats = 3
+#=========================================
+# Basic data may be set here
+number_of_measurements = 3
+number_of_repeats = 5
 t0, tf = 0.0, 10.0
 num_nodes = 500
+#=========================================
 
 m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction')
-bias_ident = sm.symbols(f'bias:{number_of_measurements}')
 
 t = me.dynamicsymbols._t
 x = sm.Matrix([me.dynamicsymbols([f'x{i}_{j}' for j in range(number_of_repeats)])
@@ -97,6 +113,8 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 
 kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction)
 
+# %%
+# Stack the equations appropriately.
 kd_total = sm.Matrix([])
 eom_total = sm.Matrix([])
 rhs_total = sm.Matrix([])
@@ -111,6 +129,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
         rhs_total = rhs_total.col_join(rhs_h)
 
 eom = kd_total.col_join(eom_total)
+
 uh =sm.Matrix([u[i, j] for i in range(number_of_measurements)
     for j in range(number_of_repeats)])
 rhs = uh.col_join(rhs_total)
@@ -123,12 +142,11 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
     sm.pprint(eom[-(number_of_repeats-i)])
 
 # %%
-# Generate Noisy Measurement Data with Bias
-# -----------------------------------------
+# Generate Noisy Measurement
+# --------------------------
 # Create 'number_of_measurements' sets of measurements with different initial
 # conditions. To get the measurements, the equations of motion are integrated,
-# and then noise is added to each point in time of the solution. Finally a bias
-# is added.
+# and then noise is added to each point in time of the solution.
 #
 
 states = [x[i, j] for i in range(number_of_measurements)
@@ -144,7 +162,8 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 
 measurements = []
 # %%
-
+# Integrate the differential equations. If np.random.seed(seed) is used, it
+# the seed must be changed for every measurement to ensure they are independent.
 for i in range(number_of_measurements):
     for j in range(number_of_repeats):
         seed = 1234*(i+1)*(j+1)
@@ -154,24 +173,10 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
                     (t0, tf), x0, t_eval=times, args=(par_vals,))
     seed = 10000 + 12345*i
     np.random.seed(seed)
-    measurements.append(sol.y[1, :] + 1.0*np.random.randn(num_nodes))
+    measurements.append(sol.y[0, :] + 1.0*np.random.randn(num_nodes))
 measurements = np.array(measurements)
-print(measurements.shape)
-print(sol.y.shape)
-
-fig, ax = plt.subplots(2*number_of_measurements*number_of_repeats, 1,
-    figsize=(6.4, 1.5*number_of_measurements*number_of_repeats),
-    sharex=True, constrained_layout=True)
-
-for i in range(2*number_of_measurements*number_of_repeats):
-    ax[i].plot(times, sol.y[i, :], color='blue')
-    ax[i].set_ylabel(f'{states[i]}')
-
-# Add bias to the measurements
-bias = 0.0 *np.random.uniform(0., 1., size=number_of_measurements)
-measurements = np.array(measurements) + bias[:, np.newaxis]
-
-print(measurements.shape)
+print('shsape of measurement array', measurements.shape)
+print('shape of solution array', sol.y.shape)
 
 # %%
 # Setup the Identification Problem
@@ -185,7 +190,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 # increased relative to the other measurements.
 #
 #
-# objective = :math:`\int_{t_0}^{t_f} \left[ \sum_{i=1}^{i = \text{number_of_measurements}} (w_i (x_i - x_i^m)^2 \right] \hspace{2pt} dt`
+# objective = :math:`\int_{t_0}^{t_f} \sum_{i=1}^{\text{number_of_measurements}} \left[ \sum_{s=1}^{\text{number_of_repeats}} (w_i (x_{i, s} - x_{i,s}^m)^2 \right] \hspace{2pt} dt`
 #
 interval_value = (tf - t0) / (num_nodes - 1)
 
@@ -197,8 +202,9 @@ def obj(free):
     sum = 0.0
     for i in range(nm):
         for j in range(nr):
-            sum += (w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] -
-                measurements[i])**2)
+            sum += np.sum((w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] -
+                measurements[i])**2))
+    return sum
 
 
 def obj_grad(free):
@@ -212,17 +218,20 @@ def obj_grad(free):
 
 
 # %%
-# By not including :math:`c`, :math:`k`, :math:`friction`, :math:`bias_i`
+# By not including :math:`c`, :math:`k`, :math:`friction`
 # in the parameter map, they will be treated as unknown parameters.
 par_map = {m: par_vals[0], l0: par_vals[3]}
 
-# %%
 bounds = {
     c: (0.01, 2.0),
     k: (0.1, 10.0),
     friction: (0.1, 10.0),
 }
-noise_scale = 2.0
+# %%
+# Set up the known trajectory map. If np.random.seed(seed) is used, the
+# seed must be changed for every map to ensure they are idependent.
+# noise_scale gives the 'strength' of the noise.
+noise_scale = 1.0
 known_trajectory_map = {}
 for i in range(number_of_measurements):
     for j in range(number_of_repeats):
@@ -230,8 +239,9 @@ def obj_grad(free):
         np.random.seed(seed)
         known_trajectory_map = (known_trajectory_map |
             {n[i, j]: noise_scale*np.random.randn(num_nodes)})
-print(len(known_trajectory_map))
-print('states', len(states))
+
+# %%
+# Set up the problem.
 problem = Problem(
     obj,
     obj_grad,
@@ -245,11 +255,9 @@ def obj_grad(free):
     bounds=bounds,
     known_trajectory_map=known_trajectory_map
 )
-print(problem.num_free)
 # %%
 # This give the sequence of the unknown parameters at the tail of solution.
 print(problem.collocator.unknown_parameters)
-print(problem.collocator.unknown_input_trajectories)
 # %%
 # Create an Initial Guess
 # -----------------------
@@ -257,7 +265,7 @@ def obj_grad(free):
 # It is reasonable to use the measurements as initial guess for the states
 # because they would be available. Here, only the measurements of the position
 # are used and the speeds are set to zero. The last values are the guesses
-# for :math:`bias_i`, :math:`c` and :math:`k`, respectively.
+# for :math:`c`, :math:`friction` and :math:`k`, respectively.
 #
 initial_guess = np.array([])
 for i in range(number_of_measurements):
@@ -267,20 +275,11 @@ def obj_grad(free):
                 np.zeros(number_of_measurements*number_of_repeats*num_nodes),
                 [0.1, 3.0, 2.0]))
 
-print(len(initial_guess))
-
 # %%
 # Solve the Optimization Problem
 # ------------------------------
 #
-# Here, the initial guess is the result of some previous run, stored in
-# 'non_contiguous_parameter_identification_bias_solution.npy'.
-#initial_guess = np.load('non_contiguous_parameter_identification_bias_solution.npy')
 solution, info = problem.solve(initial_guess)
-# %%
-# This is how the solution may be saved for a future initial guess
-# ```np.save('non_contiguous_parameter_identification_bias_solution', solution)```
-
 print(info['status_msg'])
 print(f'final value of the objective function is {info['obj_val']:.2f}' )
 
@@ -304,20 +303,35 @@ def obj_grad(free):
 
 # %%
 # Plot the measurements and the trajectories calculated.
+#-------------------------------------------------------
 #
-fig, ax = plt. subplots(number_of_measurements, 1,
-    figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True)
-for i in range(number_of_measurements):
-    ax[i].plot(times, solution[i*num_nodes:(i+1)*num_nodes], color='red')
-    ax[i].plot(times, measurements[i], color='blue', lw=0.5)
-    ax[i].set_ylabel(f'{states[i]}')
-ax[-1].set_xlabel('Time [s]')
-ax[0].set_title('Trajectories')
+sol_parsed, _, _ = parse_free(solution, len(states), 0, num_nodes )
+if number_of_measurements > 1:
+    fig, ax = plt. subplots(number_of_measurements, 1,
+        figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True)
+
+    for i in range(number_of_measurements):
+        ax[i].plot(times, sol_parsed[i*number_of_repeats, :], color='red')
+        ax[i].plot(times, measurements[i], color='blue', lw=0.5)
+        ax[i].set_ylabel(f'{states[i]}')
+    ax[-1].set_xlabel('Time [s]')
+    ax[0].set_title('Trajectories')
+
+else:
+    fig, ax = plt. subplots(1, 1, figsize=(6.4, 1.25))
+    ax.plot(times, sol_parsed[0, :], color='red')
+    ax.plot(times, measurements[0], color='blue', lw=0.5)
+    ax.set_ylabel(f'{states[0]}')
+    ax.set_xlabel('Time [s]')
+    ax.set_title('Trajectories')
+
+
 prevent_output = True
 
+# %%
+#Plot the Trajectories
+problem.plot_trajectories(solution)
 # %%
 #
 # sphinx_gallery_thumbnail_number = 3
-problem.plot_trajectories(solution)
-
 plt.show()

From 97b0da6693edecccf2754455e2c5722988efbbc6 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Wed, 25 Sep 2024 13:35:00 +0200
Subject: [PATCH 03/10] shortened title

---
 .../plot_non_contiguous_parameter_identification_wang.py  | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
index 8aeb22b7..e1cac73e 100644
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -1,6 +1,6 @@
 """
-Parameter Identification from Noncontiguous Measurements with Bias and added Noise
-==================================================================================
+Parameter Identification from Noncontiguous Measurements with Bias and Noise
+============================================================================
 
 In parameter estimation it is common to collect measurements of a system's
 trajectories from distinct experiments. For example, if you are identifying the
@@ -302,8 +302,8 @@ def obj_grad(free):
       f' Percentage error is {(solution[-2]-par_vals[-1])/solution[-2]*100:1.2f} %')
 
 # %%
-# Plot the measurements and the trajectories calculated.
-#-------------------------------------------------------
+# Plot the measurements and the trajectories calculated
+#------------------------------------------------------
 #
 sol_parsed, _, _ = parse_free(solution, len(states), 0, num_nodes )
 if number_of_measurements > 1:

From e39b75a07bcfbf269e8c28469d28ff4ed356833e Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Wed, 25 Sep 2024 13:51:34 +0200
Subject: [PATCH 04/10]  no idea which title is wrong

---
 ...t_non_contiguous_parameter_identification_wang.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
index e1cac73e..4a35ad4c 100644
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -64,13 +64,13 @@
 # %%
 # Equations of Motion.
 # --------------------
-#=========================================
+#
 # Basic data may be set here
 number_of_measurements = 3
 number_of_repeats = 5
 t0, tf = 0.0, 10.0
 num_nodes = 500
-#=========================================
+#
 
 m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction')
 
@@ -89,6 +89,7 @@
 # Form the equations of motion, such that the first number_of_repeats equations
 # belong to the first measurement, the second number_of_repeats equations belong
 # to the second measurement, and so on.
+#
 def kd_eom_rh(xh, uh, m, c, k, l0, friction):
     """" sets up the eoms for the system"""
     N = me.ReferenceFrame('N')
@@ -115,6 +116,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 
 # %%
 # Stack the equations appropriately.
+#
 kd_total = sm.Matrix([])
 eom_total = sm.Matrix([])
 rhs_total = sm.Matrix([])
@@ -164,6 +166,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 # %%
 # Integrate the differential equations. If np.random.seed(seed) is used, it
 # the seed must be changed for every measurement to ensure they are independent.
+#
 for i in range(number_of_measurements):
     for j in range(number_of_repeats):
         seed = 1234*(i+1)*(j+1)
@@ -220,6 +223,7 @@ def obj_grad(free):
 # %%
 # By not including :math:`c`, :math:`k`, :math:`friction`
 # in the parameter map, they will be treated as unknown parameters.
+#
 par_map = {m: par_vals[0], l0: par_vals[3]}
 
 bounds = {
@@ -231,6 +235,7 @@ def obj_grad(free):
 # Set up the known trajectory map. If np.random.seed(seed) is used, the
 # seed must be changed for every map to ensure they are idependent.
 # noise_scale gives the 'strength' of the noise.
+#
 noise_scale = 1.0
 known_trajectory_map = {}
 for i in range(number_of_measurements):
@@ -242,6 +247,7 @@ def obj_grad(free):
 
 # %%
 # Set up the problem.
+#
 problem = Problem(
     obj,
     obj_grad,
@@ -330,8 +336,10 @@ def obj_grad(free):
 
 # %%
 #Plot the Trajectories
+#
 problem.plot_trajectories(solution)
 # %%
 #
 # sphinx_gallery_thumbnail_number = 3
+#
 plt.show()

From f2cb12d6e038227ddc97149d9c2031fdf815e366 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Thu, 26 Sep 2024 07:59:09 +0200
Subject: [PATCH 05/10] cosmetic changes, bit faster execution

---
 ...ontiguous_parameter_identification_wang.py | 67 ++++++++++---------
 1 file changed, 35 insertions(+), 32 deletions(-)

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
index 4a35ad4c..4878433f 100644
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -1,3 +1,4 @@
+# %%
 """
 Parameter Identification from Noncontiguous Measurements with Bias and Noise
 ============================================================================
@@ -22,9 +23,9 @@
 noise corresponds to a force acting on the system).
 This idea is due to Huawei Wang and A. J. van den Bogert.
 
-So, if one has 'number_of_measurements' measurements, and one creates
-'number_of_repeats' differential equations for each measurement, one will have
-'number_of_measurements' * 'number_of_repeats' uncoupled differential equations,
+So, if one has *number_of_measurements* measurements, and one creates
+*number_of_repeats* differential equations for each measurement, one will have
+*number_of_measurements* * *number_of_repeats* uncoupled differential equations,
 all sharing the same constant parameters.
 
 Mass-spring-damper-friction Example
@@ -32,7 +33,7 @@
 
 The position of a simple system consisting of a mass connected to a fixed point
 by a spring and a damper and subject to friction is simulated and recorded as
-noisy measurements with bias. The spring constant, the damping coefficient
+noisy measurements. The spring constant, the damping coefficient
 and the coefficient of friction will be identified.
 
 **State Variables**
@@ -40,9 +41,6 @@
 - :math:`x_{i, j}`: position of the mass of the i-th measurement trial [m]
 - :math:`u_{i, j}`: speed of the mass of the i-th measurement trial [m/s]
 
-**Noise Variables**
-
-- :math:`n_{i, j}`: noise added [N]
 
 **Parameters**
 
@@ -51,6 +49,7 @@
 - :math:`k`: linear spring constant [N/m]
 - :math:`l_0`: natural length of the spring [m]
 - :math:`friction`: friction coefficient [N]
+- :math:`n_{i, j}`: noise added [N]
 
 """
 import sympy as sm
@@ -66,12 +65,12 @@
 # --------------------
 #
 # Basic data may be set here
-number_of_measurements = 3
-number_of_repeats = 5
+number_of_measurements = 4
+number_of_repeats = 3
 t0, tf = 0.0, 10.0
 num_nodes = 500
+noise_scale = 1.0
 #
-
 m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction')
 
 t = me.dynamicsymbols._t
@@ -86,9 +85,7 @@
 xh, uh = me.dynamicsymbols('xh, uh')
 
 # %%
-# Form the equations of motion, such that the first number_of_repeats equations
-# belong to the first measurement, the second number_of_repeats equations belong
-# to the second measurement, and so on.
+# Form the equations of motion, I use Kane's method here
 #
 def kd_eom_rh(xh, uh, m, c, k, l0, friction):
     """" sets up the eoms for the system"""
@@ -101,7 +98,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
     bodies = [me.Particle('part', P, m)]
     # :math:`\tanh(\alpha \cdot x) \approx sign(x)` for large :math:`\alpha`, and
     # is differentiale everywhere.
-    forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(30*uh)*N.x)]
+    forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(60*uh)*N.x)]
     kd = sm.Matrix([uh - xh.diff(t)])
     KM = me.KanesMethod(N,
                         q_ind=[xh],
@@ -115,19 +112,24 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction)
 
 # %%
-# Stack the equations appropriately.
+# Stack the equations such that the first number_of_repeats equations
+# belong to the first measurement, the second number_of_repeats equations belong
+# to the second measurement, and so on.
 #
 kd_total = sm.Matrix([])
 eom_total = sm.Matrix([])
 rhs_total = sm.Matrix([])
+
 for i in range(number_of_measurements):
     for j in range(number_of_repeats):
-        eom_h = me.msubs(fr_frstar, {xh: x[i, j], uh: u[i, j], uh.diff(t): u[i, j].diff(t)})
-        eom_h = eom_h + sm.Matrix([n[i, j]])
+        eom_h = fr_frstar.subs({xh: x[i, j], uh: u[i, j],
+            uh.diff(t): u[i, j].diff(t), xh.diff(t): u[i, j]}) + \
+            sm.Matrix([n[i, j]])
         kd_h = kd.subs({xh: x[i, j], uh: u[i, j]})
+        rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})])
+
         kd_total = kd_total.col_join(kd_h)
         eom_total = eom_total.col_join(eom_h)
-        rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})])
         rhs_total = rhs_total.col_join(rhs_h)
 
 eom = kd_total.col_join(eom_total)
@@ -155,6 +157,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
                 for j in range(number_of_repeats)] + \
         [u[i, j] for i in range(number_of_measurements)
                 for j in range(number_of_repeats)]
+
 parameters = [m, c, k, l0, friction]
 par_vals = [1.0, 0.25, 2.0, 1.0, 0.5]
 
@@ -164,22 +167,24 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction):
 
 measurements = []
 # %%
-# Integrate the differential equations. If np.random.seed(seed) is used, it
-# the seed must be changed for every measurement to ensure they are independent.
+# Integrate the differential equations. If np.random.seed(seed) is used, the
+# seed must be changed for every measurement to ensure they are independent.
 #
 for i in range(number_of_measurements):
-    for j in range(number_of_repeats):
-        seed = 1234*(i+1)*(j+1)
-        np.random.seed(seed)
-        x0 = 4*np.random.randn(2*number_of_measurements * number_of_repeats)
-        sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(),
-                    (t0, tf), x0, t_eval=times, args=(par_vals,))
-    seed = 10000 + 12345*i
+    seed = 1234*(i+1)
+    np.random.seed(seed)
+
+    x0 = 3.0*np.random.randn(2*number_of_measurements * number_of_repeats)
+    sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(),
+            (t0, tf), x0, t_eval=times, args=(par_vals,))
+
+    seed = 10000 + 12345*(i+1)
     np.random.seed(seed)
     measurements.append(sol.y[0, :] + 1.0*np.random.randn(num_nodes))
+
 measurements = np.array(measurements)
-print('shsape of measurement array', measurements.shape)
-print('shape of solution array', sol.y.shape)
+print('shape of measurement array', measurements.shape)
+print('shape of solution array   ', sol.y.shape)
 
 # %%
 # Setup the Identification Problem
@@ -209,7 +214,6 @@ def obj(free):
                 measurements[i])**2))
     return sum
 
-
 def obj_grad(free):
     grad = np.zeros_like(free)
     for i in range(nm):
@@ -236,7 +240,6 @@ def obj_grad(free):
 # seed must be changed for every map to ensure they are idependent.
 # noise_scale gives the 'strength' of the noise.
 #
-noise_scale = 1.0
 known_trajectory_map = {}
 for i in range(number_of_measurements):
     for j in range(number_of_repeats):
@@ -319,7 +322,7 @@ def obj_grad(free):
     for i in range(number_of_measurements):
         ax[i].plot(times, sol_parsed[i*number_of_repeats, :], color='red')
         ax[i].plot(times, measurements[i], color='blue', lw=0.5)
-        ax[i].set_ylabel(f'{states[i]}')
+        ax[i].set_ylabel(f'{states[number_of_repeats* i]}')
     ax[-1].set_xlabel('Time [s]')
     ax[0].set_title('Trajectories')
 

From 3523f76afe1cdad6476a8a2cdb9c5e0be43c730d Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Thu, 26 Sep 2024 08:08:52 +0200
Subject: [PATCH 06/10] try again

---
 .../plot_non_contiguous_parameter_identification_wang.py        | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
index 4878433f..2597502f 100644
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -85,7 +85,7 @@
 xh, uh = me.dynamicsymbols('xh, uh')
 
 # %%
-# Form the equations of motion, I use Kane's method here
+# Form the equations of motion, I use Kane's method here.
 #
 def kd_eom_rh(xh, uh, m, c, k, l0, friction):
     """" sets up the eoms for the system"""

From c0cd80599d81bb3c7ee0faa31c28bf6b9b0855c3 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Sat, 16 Nov 2024 14:44:38 +0100
Subject: [PATCH 07/10] just to see if it wasses all test now

---
 .../plot_non_contiguous_parameter_identification_wang.py         | 1 -
 1 file changed, 1 deletion(-)

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
index 2597502f..fd96d8c9 100644
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
@@ -1,4 +1,3 @@
-# %%
 """
 Parameter Identification from Noncontiguous Measurements with Bias and Noise
 ============================================================================

From 3f6139d5bf02f62dd0637e47b856d7b1cdb62320 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Tue, 19 Nov 2024 10:49:12 +0100
Subject: [PATCH 08/10] betts examples 10.144 and 10.145

---
 examples-gallery/plot_betts_10_144_145.py | 220 ++++++++++++++++++++++
 1 file changed, 220 insertions(+)
 create mode 100644 examples-gallery/plot_betts_10_144_145.py

diff --git a/examples-gallery/plot_betts_10_144_145.py b/examples-gallery/plot_betts_10_144_145.py
new file mode 100644
index 00000000..176bdbe3
--- /dev/null
+++ b/examples-gallery/plot_betts_10_144_145.py
@@ -0,0 +1,220 @@
+"""
+Van der Pol Oscillator
+======================
+
+https://en.wikipedia.org/wiki/Van_der_Pol_oscillator
+
+These are example 10.144 / 145 from Betts' book "Practical Methods for Optimal Control
+Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems.
+It is described in more detail in section 4.14. example 4.11 of the book.
+
+"""
+import numpy as np
+import sympy as sm
+import sympy.physics.mechanics as me
+from opty.direct_collocation import Problem
+from opty.utils import create_objective_function
+import matplotlib.pyplot as plt
+
+# %%
+# First version of the problem (10.144).
+# --------------------------------------
+# **States**
+# :math:`y_1, y_2` : state variables
+#
+# **Controls**
+# :math:`u` : control variables
+#
+# Equations of Motion.
+# --------------------
+t = me.dynamicsymbols._t
+
+y1, y2 = me.dynamicsymbols('y1, y2')
+u = me.dynamicsymbols('u')
+
+eom = sm.Matrix([
+    -y1.diff() + y2,
+    -y2.diff(t) + (1 - y1**2)*y2 - y1 + u,
+])
+sm.pprint(eom)
+
+# %%
+# Define and Solve the Optimization Problem.
+# -----------------------------------------
+num_nodes = 2001
+iterations = 1000
+
+t0, tf = 0.0, 5.0
+interval_value = (tf - t0)/(num_nodes - 1)
+
+state_symbols = (y1, y2)
+unkonwn_input_trajectories = (u, )
+
+objective = sm.integrate(u**2 + y1**2 + y2**2, t)
+obj, obj_grad = create_objective_function(objective,
+                state_symbols,
+                unkonwn_input_trajectories,
+                tuple(),
+                num_nodes,
+                interval_value
+)
+
+instance_constraints = (
+        y1.func(t0) - 1.0,
+        y2.func(t0),
+)
+
+bounds = {
+        y2: (-0.4, np.inf),
+}
+
+prob = Problem(
+        obj,
+        obj_grad,
+        eom,
+        state_symbols,
+        num_nodes,
+        interval_value,
+        instance_constraints= instance_constraints,
+        bounds=bounds,
+)
+
+# %%
+# Solve the optimization problem.
+# Give some rough estimates for the trajectories.
+
+initial_guess = np.ones(prob.num_free)
+
+# Find the optimal solution.
+for i in range(1):
+    solution, info = prob.solve(initial_guess)
+    initial_guess = solution
+    print(info['status_msg'])
+    print(f'Objectve is: {info['obj_val']:.8f}, ' +
+          f'as per the book it is {2.95369916}, so the deviation is: ' +
+        f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %')
+solution1 = solution
+# %%
+# Plot the optimal state and input trajectories.
+prob.plot_trajectories(solution)
+
+# %%
+# Plot the constraint violations.
+prob.plot_constraint_violations(solution)
+
+# %%
+# Plot the objective function as a function of optimizer iteration.
+prob.plot_objective_value()
+
+# %%
+# Second version of the problem (10.145).
+# ---------------------------------------
+# This is example 10.145 from Betts' book "Practical Methods for Optimal Control
+# Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems.
+# It is same as problem 10.144 but a bit reformulated.
+# It has two control variables and one additional algebraic equation of motion.
+# As opty needs as many state variables as equations of motion, I simply call
+# the additional control variable v a state variable.
+# It is described in more detail in section 4.14, example 4.11 of the book.
+#
+# This formulation seems to be more accurate compared to the above, when
+# considering the violations of the constraints.
+#
+# **States**
+# :math:`y_1, y_2, v` : state variables
+#
+# **Controls**
+# :math:`u` : control variables
+#
+# Equations of Motion.
+# --------------------
+y1, y2, v = me.dynamicsymbols('y1, y2, v')
+u = me.dynamicsymbols('u')
+
+eom = sm.Matrix([
+    -y1.diff() + y2,
+    -y2.diff(t) + v - y1 + u,
+    v - (1-y1**2)*y2,
+])
+sm.pprint(eom)
+
+# %%
+# Define and Solve the Optimization Problem.
+# ------------------------------------------
+state_symbols = (y1, y2, v)
+unkonwn_input_trajectories = (u,)
+
+objective = sm.integrate(u**2 + y1**2 + y2**2, t)
+obj, obj_grad = create_objective_function(objective,
+                state_symbols,
+                unkonwn_input_trajectories,
+                tuple(),
+                num_nodes,
+                interval_value
+)
+
+instance_constraints = (
+        y1.func(t0) - 1.0,
+        y2.func(t0),
+)
+
+bounds = {
+        y2: (-0.4, np.inf),
+}
+
+prob = Problem(
+        obj,
+        obj_grad,
+        eom,
+        state_symbols,
+        num_nodes,
+        interval_value,
+        instance_constraints= instance_constraints,
+        bounds=bounds,
+)
+
+# %%
+# Solve the optimization problem.
+# Give some rough estimates for the trajectories.
+
+initial_guess = np.ones(prob.num_free)
+
+# Find the optimal solution.
+for i in range(1):
+    solution, info = prob.solve(initial_guess)
+    initial_guess = solution
+    print(info['status_msg'])
+    print(f'Objectve is: {info['obj_val']:.8f}, ' +
+          f'as per the book it is {2.95369916}, so the deviation is: ' +
+        f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %')
+# %%
+# Plot the optimal state and input trajectories.
+prob.plot_trajectories(solution)
+
+# %%
+# Plot the constraint violations.
+prob.plot_constraint_violations(solution)
+
+# %%
+# Plot the objective function as a function of optimizer iteration.
+prob.plot_objective_value()
+
+# %%
+# Plot the Difference between the two Solutions.
+# ----------------------------------------------
+diffy1 = solution1[: num_nodes] - solution[: num_nodes]
+diffy2 = solution1[num_nodes: 2*num_nodes] - solution[num_nodes: 2*num_nodes]
+diffu = solution1[2*num_nodes:] - solution[3*num_nodes:]
+times = np.linspace(t0, tf, num_nodes)
+
+fig, ax = plt.subplots(3, 1, figsize=(7, 4), sharex=True, constrained_layout=True)
+ax[0].plot(times, diffy1, label='Delta y1')
+ax[0].legend()
+ax[1].plot(times, diffy2, label='Delta y2')
+ax[1].legend()
+ax[2].plot(times, diffu, label='Delta u')
+ax[2].legend()
+ax[2].set_xlabel('Time')
+ax[0].set_title('Differences between the two solutions')
+avoid_printing = True
+# sphinx_gallery_thumbnail_number = 2

From 4c60b6794ba9c72e122a7dfba0103fd6faf59701 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <83544146+Peter230655@users.noreply.github.com>
Date: Tue, 19 Nov 2024 11:02:36 +0100
Subject: [PATCH 09/10] Delete
 examples-gallery/plot_non_contiguous_parameter_identification_wang.py

---
 ...ontiguous_parameter_identification_wang.py | 347 ------------------
 1 file changed, 347 deletions(-)
 delete mode 100644 examples-gallery/plot_non_contiguous_parameter_identification_wang.py

diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
deleted file mode 100644
index fd96d8c9..00000000
--- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py
+++ /dev/null
@@ -1,347 +0,0 @@
-"""
-Parameter Identification from Noncontiguous Measurements with Bias and Noise
-============================================================================
-
-In parameter estimation it is common to collect measurements of a system's
-trajectories from distinct experiments. For example, if you are identifying the
-parameters of a mass-spring-damper system, you may excite the system with
-different initial conditions multiple times and measure the position of the
-mass. The data cannot simply be stacked end-to-end in time and the
-identification run because the measurement data would be discontinuous between
-each measurement trial.
-
-A workaround in opty is to create a set of differential equations with unique
-state variables for each measurement trial that all share the same constant
-parameters. You can then identify the parameters from all measurement trials
-simultaneously by passing the uncoupled differential equations to opty.
-This idea is due to Jason Moore.
-
-In addition, for each measurement, one may create a set of differential
-equations with unique state variables that again share the same constant
-parameters. In addition, one adds noise to these equations. (Intuitively this
-noise corresponds to a force acting on the system).
-This idea is due to Huawei Wang and A. J. van den Bogert.
-
-So, if one has *number_of_measurements* measurements, and one creates
-*number_of_repeats* differential equations for each measurement, one will have
-*number_of_measurements* * *number_of_repeats* uncoupled differential equations,
-all sharing the same constant parameters.
-
-Mass-spring-damper-friction Example
-===================================
-
-The position of a simple system consisting of a mass connected to a fixed point
-by a spring and a damper and subject to friction is simulated and recorded as
-noisy measurements. The spring constant, the damping coefficient
-and the coefficient of friction will be identified.
-
-**State Variables**
-
-- :math:`x_{i, j}`: position of the mass of the i-th measurement trial [m]
-- :math:`u_{i, j}`: speed of the mass of the i-th measurement trial [m/s]
-
-
-**Parameters**
-
-- :math:`m`: mass [kg]
-- :math:`c`: linear damping coefficient [Ns/m]
-- :math:`k`: linear spring constant [N/m]
-- :math:`l_0`: natural length of the spring [m]
-- :math:`friction`: friction coefficient [N]
-- :math:`n_{i, j}`: noise added [N]
-
-"""
-import sympy as sm
-import sympy.physics.mechanics as me
-import numpy as np
-from scipy.integrate import solve_ivp
-from opty import Problem
-from opty.utils import parse_free
-import matplotlib.pyplot as plt
-
-# %%
-# Equations of Motion.
-# --------------------
-#
-# Basic data may be set here
-number_of_measurements = 4
-number_of_repeats = 3
-t0, tf = 0.0, 10.0
-num_nodes = 500
-noise_scale = 1.0
-#
-m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction')
-
-t = me.dynamicsymbols._t
-x = sm.Matrix([me.dynamicsymbols([f'x{i}_{j}' for j in range(number_of_repeats)])
-        for i in range(number_of_measurements)])
-u = sm.Matrix([me.dynamicsymbols([f'u{i}_{j}' for j in range(number_of_repeats)])
-        for i in range(number_of_measurements)])
-
-n = sm.Matrix([me.dynamicsymbols([f'n{i}_{j}' for j in range(number_of_repeats)])
-        for i in range(number_of_measurements)])
-
-xh, uh = me.dynamicsymbols('xh, uh')
-
-# %%
-# Form the equations of motion, I use Kane's method here.
-#
-def kd_eom_rh(xh, uh, m, c, k, l0, friction):
-    """" sets up the eoms for the system"""
-    N = me.ReferenceFrame('N')
-    O, P = sm.symbols('O, P', cls=me.Point)
-    O.set_vel(N, 0)
-
-    P.set_pos(O, xh*N.x)
-    P.set_vel(N, uh*N.x)
-    bodies = [me.Particle('part', P, m)]
-    # :math:`\tanh(\alpha \cdot x) \approx sign(x)` for large :math:`\alpha`, and
-    # is differentiale everywhere.
-    forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(60*uh)*N.x)]
-    kd = sm.Matrix([uh - xh.diff(t)])
-    KM = me.KanesMethod(N,
-                        q_ind=[xh],
-                        u_ind=[uh],
-                        kd_eqs=kd
-    )
-    fr, frstar = KM.kanes_equations(bodies, forces)
-    rhs = KM.rhs()
-    return (kd, fr + frstar, rhs)
-
-kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction)
-
-# %%
-# Stack the equations such that the first number_of_repeats equations
-# belong to the first measurement, the second number_of_repeats equations belong
-# to the second measurement, and so on.
-#
-kd_total = sm.Matrix([])
-eom_total = sm.Matrix([])
-rhs_total = sm.Matrix([])
-
-for i in range(number_of_measurements):
-    for j in range(number_of_repeats):
-        eom_h = fr_frstar.subs({xh: x[i, j], uh: u[i, j],
-            uh.diff(t): u[i, j].diff(t), xh.diff(t): u[i, j]}) + \
-            sm.Matrix([n[i, j]])
-        kd_h = kd.subs({xh: x[i, j], uh: u[i, j]})
-        rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})])
-
-        kd_total = kd_total.col_join(kd_h)
-        eom_total = eom_total.col_join(eom_h)
-        rhs_total = rhs_total.col_join(rhs_h)
-
-eom = kd_total.col_join(eom_total)
-
-uh =sm.Matrix([u[i, j] for i in range(number_of_measurements)
-    for j in range(number_of_repeats)])
-rhs = uh.col_join(rhs_total)
-
-for i in range(number_of_repeats):
-    sm.pprint(eom[i])
-for i in range(3):
-    print('.........')
-for i in range(number_of_repeats):
-    sm.pprint(eom[-(number_of_repeats-i)])
-
-# %%
-# Generate Noisy Measurement
-# --------------------------
-# Create 'number_of_measurements' sets of measurements with different initial
-# conditions. To get the measurements, the equations of motion are integrated,
-# and then noise is added to each point in time of the solution.
-#
-
-states = [x[i, j] for i in range(number_of_measurements)
-                for j in range(number_of_repeats)] + \
-        [u[i, j] for i in range(number_of_measurements)
-                for j in range(number_of_repeats)]
-
-parameters = [m, c, k, l0, friction]
-par_vals = [1.0, 0.25, 2.0, 1.0, 0.5]
-
-eval_rhs = sm.lambdify(states + parameters, rhs)
-
-times = np.linspace(t0, tf, num=num_nodes)
-
-measurements = []
-# %%
-# Integrate the differential equations. If np.random.seed(seed) is used, the
-# seed must be changed for every measurement to ensure they are independent.
-#
-for i in range(number_of_measurements):
-    seed = 1234*(i+1)
-    np.random.seed(seed)
-
-    x0 = 3.0*np.random.randn(2*number_of_measurements * number_of_repeats)
-    sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(),
-            (t0, tf), x0, t_eval=times, args=(par_vals,))
-
-    seed = 10000 + 12345*(i+1)
-    np.random.seed(seed)
-    measurements.append(sol.y[0, :] + 1.0*np.random.randn(num_nodes))
-
-measurements = np.array(measurements)
-print('shape of measurement array', measurements.shape)
-print('shape of solution array   ', sol.y.shape)
-
-# %%
-# Setup the Identification Problem
-# --------------------------------
-#
-# The goal is to identify the damping coefficient :math:`c`, the spring
-# constant :math:`k` and the coefficient of friction :math:`friction`.
-# The objective :math:`J` is to minimize the least square
-# difference in the optimal simulation as compared to the measurements.  If
-# some measurement is considered more reliable, its weight :math:`w` may be
-# increased relative to the other measurements.
-#
-#
-# objective = :math:`\int_{t_0}^{t_f} \sum_{i=1}^{\text{number_of_measurements}} \left[ \sum_{s=1}^{\text{number_of_repeats}} (w_i (x_{i, s} - x_{i,s}^m)^2 \right] \hspace{2pt} dt`
-#
-interval_value = (tf - t0) / (num_nodes - 1)
-
-w =[1 for _ in range(number_of_measurements)]
-
-nm = number_of_measurements
-nr = number_of_repeats
-def obj(free):
-    sum = 0.0
-    for i in range(nm):
-        for j in range(nr):
-            sum += np.sum((w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] -
-                measurements[i])**2))
-    return sum
-
-def obj_grad(free):
-    grad = np.zeros_like(free)
-    for i in range(nm):
-        for j in range(nr):
-            grad[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] = 2*w[i]*interval_value*(
-                free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - measurements[i]
-            )
-    return grad
-
-
-# %%
-# By not including :math:`c`, :math:`k`, :math:`friction`
-# in the parameter map, they will be treated as unknown parameters.
-#
-par_map = {m: par_vals[0], l0: par_vals[3]}
-
-bounds = {
-    c: (0.01, 2.0),
-    k: (0.1, 10.0),
-    friction: (0.1, 10.0),
-}
-# %%
-# Set up the known trajectory map. If np.random.seed(seed) is used, the
-# seed must be changed for every map to ensure they are idependent.
-# noise_scale gives the 'strength' of the noise.
-#
-known_trajectory_map = {}
-for i in range(number_of_measurements):
-    for j in range(number_of_repeats):
-        seed = 10000 + 12345*(i+1)*(j+1)
-        np.random.seed(seed)
-        known_trajectory_map = (known_trajectory_map |
-            {n[i, j]: noise_scale*np.random.randn(num_nodes)})
-
-# %%
-# Set up the problem.
-#
-problem = Problem(
-    obj,
-    obj_grad,
-    eom,
-    states,
-    num_nodes,
-    interval_value,
-    known_parameter_map=par_map,
-    time_symbol=me.dynamicsymbols._t,
-    integration_method='midpoint',
-    bounds=bounds,
-    known_trajectory_map=known_trajectory_map
-)
-# %%
-# This give the sequence of the unknown parameters at the tail of solution.
-print(problem.collocator.unknown_parameters)
-# %%
-# Create an Initial Guess
-# -----------------------
-#
-# It is reasonable to use the measurements as initial guess for the states
-# because they would be available. Here, only the measurements of the position
-# are used and the speeds are set to zero. The last values are the guesses
-# for :math:`c`, :math:`friction` and :math:`k`, respectively.
-#
-initial_guess = np.array([])
-for i in range(number_of_measurements):
-    for j in range(number_of_repeats):
-        initial_guess = np.concatenate((initial_guess, measurements[i, :]))
-initial_guess = np.hstack((initial_guess,
-                np.zeros(number_of_measurements*number_of_repeats*num_nodes),
-                [0.1, 3.0, 2.0]))
-
-# %%
-# Solve the Optimization Problem
-# ------------------------------
-#
-solution, info = problem.solve(initial_guess)
-print(info['status_msg'])
-print(f'final value of the objective function is {info['obj_val']:.2f}' )
-
-# %%
-problem.plot_objective_value()
-
-# %%
-problem.plot_constraint_violations(solution)
-
-
-# %%
-# The identified parameters are:
-# ------------------------------
-#
-print(f'Estimate of damping coefficient is      {solution[-3]: 1.2f}' +
-      f' Percentage error is {(solution[-3]-par_vals[1])/solution[-3]*100:1.2f} %')
-print(f'Estimate of the spring constant is      {solution[-1]: 1.2f}' +
-      f' Percentage error is {(solution[-1]-par_vals[2])/solution[-1]*100:1.2f} %')
-print(f'Estimate of the friction coefficient is {solution[-2]: 1.2f}'
-      f' Percentage error is {(solution[-2]-par_vals[-1])/solution[-2]*100:1.2f} %')
-
-# %%
-# Plot the measurements and the trajectories calculated
-#------------------------------------------------------
-#
-sol_parsed, _, _ = parse_free(solution, len(states), 0, num_nodes )
-if number_of_measurements > 1:
-    fig, ax = plt. subplots(number_of_measurements, 1,
-        figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True)
-
-    for i in range(number_of_measurements):
-        ax[i].plot(times, sol_parsed[i*number_of_repeats, :], color='red')
-        ax[i].plot(times, measurements[i], color='blue', lw=0.5)
-        ax[i].set_ylabel(f'{states[number_of_repeats* i]}')
-    ax[-1].set_xlabel('Time [s]')
-    ax[0].set_title('Trajectories')
-
-else:
-    fig, ax = plt. subplots(1, 1, figsize=(6.4, 1.25))
-    ax.plot(times, sol_parsed[0, :], color='red')
-    ax.plot(times, measurements[0], color='blue', lw=0.5)
-    ax.set_ylabel(f'{states[0]}')
-    ax.set_xlabel('Time [s]')
-    ax.set_title('Trajectories')
-
-
-prevent_output = True
-
-# %%
-#Plot the Trajectories
-#
-problem.plot_trajectories(solution)
-# %%
-#
-# sphinx_gallery_thumbnail_number = 3
-#
-plt.show()

From 936aebe6a887bdd5c4d4add11f5f935f88d7d229 Mon Sep 17 00:00:00 2001
From: Peter Stahlecker <peter.stahlecker@gmail.com>
Date: Tue, 19 Nov 2024 11:42:27 +0100
Subject: [PATCH 10/10] changed underline of a title

---
 examples-gallery/plot_betts_10_144_145.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples-gallery/plot_betts_10_144_145.py b/examples-gallery/plot_betts_10_144_145.py
index 176bdbe3..6c626b63 100644
--- a/examples-gallery/plot_betts_10_144_145.py
+++ b/examples-gallery/plot_betts_10_144_145.py
@@ -40,7 +40,7 @@
 
 # %%
 # Define and Solve the Optimization Problem.
-# -----------------------------------------
+# ------------------------------------------
 num_nodes = 2001
 iterations = 1000