# Introduction

Given a closed loop system and a simpler identification model of the system we want to estimate the parameters in the controller block of the closed loop system.

Here I try to show that a systematic error in the identified parameters exists and that is a function of the ratio of the external pertubation and the process noise. This bias is present in indirect identification and direct identification methods.

Assumptions:
- Both the true system and the identification model have the same states.
- We can model the plant perfectly in the identification model.

This import makes some of the helper functions in the `src` directory available.


```
import sys
sys.path.append('../src')
```

# Derivation of the Models

The following cell derives the closed-loop non-linear equations of motion for a point mass pendulum on a laterally moving base, i.e. the classic inverted pendulum on a moving cart, which controlled by a simple proportional torque controller at the pin joint and subjected to reference noise, $w$.

```
  iy
  ^       E
  |    |   o m
  |    |0 /      |
  |    | / l     v g
  |   J|/
O -----o------> ix
  |-x->|
```

Symbols:
- $g$: acceleration due to gravity
- $m$: bob mass
- $l$: pendulum length
- $k$: proportional controller gain
- $\hat{\mathbf{i}}_{xyz}$: inertial frame unit vectors
- $\hat{\mathbf{p}}_{xyz}$: pendulum frame unit vectors
- $O$: point that locates the pin joint
- $B$: point that locates the bob
- $x(t), v(t), a(t)$: base position, velocity, and acceleration
- $\theta(t), \omega(t)$: pendulum angle and angular rate
- $w(t)$: white noise on the angle reference tracking


```
import sympy as sm
import sympy.physics.mechanics as me
```


```
sm.physics.vector.init_vprinting()
```


```
m, l, g, k = sm.symbols('m, l, g, k')
theta, omega, x, v, a, w = me.dynamicsymbols('theta, omega, x, v, a, w')

inertial_frame = me.ReferenceFrame('I')
pendulum_frame = inertial_frame.orientnew('P', 'Axis', (theta,
                                                        inertial_frame.z))

pendulum_frame.set_ang_vel(inertial_frame, omega * inertial_frame.z)

origin = me.Point('O')

joint = origin.locatenew('J', x * inertial_frame.x)
joint.set_vel(inertial_frame, v * inertial_frame.x)
joint.set_acc(inertial_frame, a * inertial_frame.x)

end = joint.locatenew('E', l * pendulum_frame.y)
end.v2pt_theory(joint, inertial_frame, pendulum_frame)
end.a2pt_theory(joint, inertial_frame, pendulum_frame)

bob = me.Particle('B', end, m)

kane = me.KanesMethod(inertial_frame, (theta,), (omega,),
                      (theta.diff() - omega,))

gravitational_force = -m * g * inertial_frame.y

control_torque = k * (-theta + w) * inertial_frame.z

loads = [(end, gravitational_force), (pendulum_frame, control_torque)]

fr, frstar = kane.kanes_equations(loads, [bob])
```


```
x = kane.q.col_join(kane.u)
xdot = x.diff()
```

## Non-linear Model

The explicit first order non-linear ordinary differential equations take the following form:


```
non_linear_sym_rhs = sm.simplify(kane.rhs())
sm.Eq(xdot, non_linear_sym_rhs.expand())
```

The model above is what I'll call the true system. We will use a slightly different model to identify parameters of this true system. The identification model will simply be the above model with $w=0$ and I'll introduce new variable names for the parameters that we will identify, in this case simply the gain.


```
k_id = sm.symbols('k_id')
non_linear_sym_id_rhs = non_linear_sym_rhs.subs({w: 0, k: k_id}).expand()
sm.Eq(xdot, non_linear_sym_id_rhs)
```

Now, if $k_{id}$ is the only free parameter in the identifaction process and we use data, with no measurement noise, generated from integrating the true system, then $k_{id}$ will be adjusted such that it accounts for the motion due to the noise term in the model: $\frac{kw}{l^2m}$. By equating the two models and solving for $k_{id}$ we can how the true gain relates to the gain we identify.


```
sol = sm.solve(non_linear_sym_rhs - non_linear_sym_id_rhs, k_id)
sm.Eq(k_id, sm.factor_terms(sol[k_id]))
```

## Linear Model

The system with a linear plant gives the same result:


```
A, B, u = kane.linearize(op_point={theta: 0, omega: 0}, new_method=True, A_and_B=True)
```


```
linear_sym_rhs = sm.simplify(A * x + B * u)
sm.Eq(xdot, linear_sym_rhs.expand())
```


```
linear_sym_id_rhs = linear_sym_rhs.subs({w: 0, k: k_id}).expand()
sm.Eq(xdot, linear_sym_id_rhs)
```


```
sol = sm.solve(linear_sym_rhs - linear_sym_id_rhs, k_id)
sm.Eq(k_id, sm.factor_terms(sol[k_id]))
```

So this tells us that for $k_{id}$ to equal $k$, $\frac{w}{\theta} \rightarrow 0$. So, $w << \theta$ and that can happen if $w$ is very small relative to $\theta$. If $w$ is small it is equivalent to saying that our identification model is a good model of the process. But if it is unknown how bad our model is, $w$ has some magnitude, i.e. the modeling error. But this equation says that if we can cause $\theta$ to be significantly larger than $w$, then $k_{id} \rightarrow k$. $w$ and $\theta$ are both time varying parameters, so it is difficult to say what $\frac{w}{\theta}$ is and they could certainly be correlated. For example, if we choose a linearized model as our identification model, then surely the process noise, $w$, is correlated to $\theta$, i.e. if you exercise the linear model in the non-linear regime, i.e. large $\theta$, then the process noise term must account for this $\theta$ dependent error. But, if they aren't necessarily correlated, one way to increase $\theta$ without increasing $w$ in this model is to increase $a$. By perturbing the model with a term that isn't subjected to the feedback loop, we can also potentially have $k_id \rightarrow k$.

So, the conclusion here is that the ratio of $\frac{a}{w}$ must be large for $k_id \rightarrow k$. This is an example of systematic error due to our identification model being a poor estimate of a true process. I'm not sure how to generalize this yet, but I'm pretty sure it is true for other systems.

Ljung 1999 Chapter 13 sections 4-5 (not so clearly) state that there are bias issues when doing indirect identification, just like there are with direct indetification. The bias is due to a number of things, but one is that you don't have a good noise model, i.e. that your identification model doesn't represent the true system well. There are ways to reduce the bias, of course, but for complex systems where little knowledge is known about any of the processes (i.e. the human) it is difficult to ensure that you are reducing this bias.

Some notes from Ljung 1999 that are relevant:

- pg 432: "The closed loop experiment is informative if and only if $r(t)$ is persistently excited." $r(t)$ is the reference signal in a standard closed loop system when you want to identify the plant. So I think, for our case, this means that the platform acceleration must be persistently excited.
- The previous statement is for linear systems, so he clarifies that "nonlinear of time varying or complex (high order) regulators should in general yield experiments that are informative enough." The regulator is the plant in our case, so we have a complex plant (the human open loop system). So if you can't excite the platform enough, you may still have informative data due to the non-linearity of the plant.
- pg 433: "the bias inclination will be small in freqeuncy ranges where either (or all) of the following holds: 1) the noise model is good, 2) the feedback contribution to the input spectrum is small, 3) the signal to noise ratio is good." In terms of our problem 1) refers to the process noise, i.e. if you don't include a good noise model in your identification model, you get biased results. This is what I show above. I think 2 and 3 essentially mean, in our case, that the reference noise is small relative to the input to the controlle.
- pg 437 shows that you get the same bias as the direct approach (13.62)

# Simulate the True System


```
from collections import OrderedDict

import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from pydy.codegen.ode_function_generators import generate_ode_function

from fast_interpolate import Interpolator
```

First I'll pick some numerical values for the constant parameters.


```
par_map = OrderedDict()
par_map[m] = 1.0
par_map[l] = 1.0
par_map[g] = 9.81
par_map[k] = 15.0
par_map
```

Here I generate numerical functions that can evaluate the right hand side of the closed loop ordinary differential equations for both the non-linear and linear models.


```
non_linear_rhs = generate_ode_function(non_linear_sym_rhs,
                                       kane.q, kane.u, par_map.keys(),
                                       specifieds=u,
                                       constants_arg_type='array',
                                       specifieds_arg_type='function',
                                       generator='cython')
```


```
linear_rhs = generate_ode_function(linear_sym_rhs,
                                   kane.q, kane.u, par_map.keys(),
                                   specifieds=u,
                                   constants_arg_type='array',
                                   specifieds_arg_type='function',
                                   generator='cython')
```

Now I pick a simulation duration and the number of time points.


```
num_nodes = 1001
duration = 10.0

time = np.linspace(0.0, duration, num_nodes)

interval = duration / (num_nodes - 1)
```

Next, I specify the time histories of the inputs so that we can simulate the system. I'll use a sum of sines for both the acceleration $a(t)$ and the reference noise $w(t)$.


```
def sum_of_sines(time, amplitudes, frequencies):
    
    phases = 2.0 * np.pi * np.random.ranf(len(frequencies))

    sines = np.zeros_like(time)

    for a, w, p in zip(amplitudes, frequencies, phases):
        sines += a * np.sin(w * time + p)

    return sines
```


```
u
```


```
frequencies = np.logspace(np.log10(0.15), np.log10(15.0), num=12)

specifieds = np.empty((len(time), 2), dtype=float)

acceleration = sum_of_sines(time, 0.2 * np.ones_like(frequencies), frequencies)
process_noise = sum_of_sines(time, 0.005 * np.ones_like(frequencies), frequencies)

specifieds[:, 0] = acceleration
specifieds[:, 1] = process_noise
```

Because I will use a variable step integrator, I create a function that interpolates the specified inputs at any time.


```
interpolator = Interpolator(time, specifieds)

result = np.zeros(2)

def perturb(x, t):
    if t > time[-1]:
        return interpolator.interpolate(time[-1], result)
    else:
        return interpolator.interpolate(t, result)
```

Initial state is zero. I make it a small number here to avoid a divide by zero later on.


```
x0 = 0.0001 * np.ones(2)
```

Now I simulate the non-linear and linear systems under the influence of the two specified inputs.


```
x_non_lin_traj = odeint(non_linear_rhs, x0, time,
                        args=(perturb, np.array(par_map.values())))
```


```
x_lin_traj = odeint(linear_rhs, x0, time,
                    args=(perturb, np.array(par_map.values())))
```

And plot the resulting state trajectories.


```
%matplotlib inline
```


```
fig, axes = plt.subplots(4, 1, sharex=True, figsize=(8, 8))

axes[0].plot(time, acceleration)
axes[0].set_ylabel(r'a(t) [$ms^{-2}$]')

axes[1].plot(time, np.rad2deg(process_noise))
axes[1].set_ylabel(r'w(t) [$deg$]')

axes[2].plot(time, np.rad2deg(x_non_lin_traj[:, 0]), 'red', label='non-linear')
axes[2].plot(time, np.rad2deg(x_lin_traj[:, 0]), 'blue', label='linear')
axes[2].set_ylabel(r'$\theta(t)$ [deg]')

axes[3].plot(time, np.rad2deg(x_non_lin_traj[:, 1]), 'red', label='non-linear')
axes[3].plot(time, np.rad2deg(x_lin_traj[:, 1]), 'blue', label='linear')
axes[3].set_ylabel(r'$\omega(t)$ [deg/s]')
axes[3].set_xlabel('Time [s]')

leg = plt.legend()
```

## The Ratio

I'm curious what $k (1 - \frac{w}{\theta})$ looks like. It should hover around $k$ if we expect the identified gain to be anywhere close.


```
plt.plot(time, par_map[k] * (1 - process_noise / x_non_lin_traj[:, 0]))
plt.ylim([par_map[k] - 5.0, par_map[k] + 5.0])
```

The graph doesn't give a good impression but if you divided the standard devations of the process noise and $\theta$ I get a number that may tell something useful:


```
par_map[k] * (1 - process_noise.std() / x_non_lin_traj[:, 0].std())
```

# Indirect Parameter Identification with Direct Collocation


```
from opty.direct_collocation import Problem
```

Now I setup a direct collocation problem that simply minimizes the difference in the identification model's trajectory and the simulated trajectories from above. I'll start by using the non-linear identification model.


```
def obj(free):
    return interval * np.sum((x_non_lin_traj.T.flatten() - free[:-1])**2)

def obj_grad(free):
    grad = np.zeros_like(free)
    grad[:-1] = 2.0 * interval * (free[:-1] - x_non_lin_traj.T.flatten())
    return grad
```


```
known_par_map = par_map.copy()
known_par_map.pop(k)
```

## Non-linear Model


```
prob = Problem(obj,
               obj_grad,
               xdot - non_linear_sym_id_rhs,
               x,
               num_nodes,
               interval,
               known_parameter_map=known_par_map,
               known_trajectory_map={a: acceleration},
               integration_method='midpoint')
```


```
prob.addOption('linear_solver', 'ma57')
```

I'll give the best initial guess I can, just for quick and easy convergence. The "measured" states and the known gain are given.


```
initial_guess = np.hstack((x_non_lin_traj.T.flatten(), par_map[k]))
```

And the solution:


```
sol, info = prob.solve(initial_guess)
sm.Eq(k_id, sol[-1])
```

This compares the actual trajectory to the direct collocation trajectory solution.


```
fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8, 8))

axes[0].plot(time, acceleration)
axes[0].set_ylabel(r'a(t) [$ms^{-2}$]')

axes[1].plot(time, np.rad2deg(x_non_lin_traj[:, 0]), '.k')
axes[1].plot(time, np.rad2deg(sol[:-1].reshape(2, num_nodes).T)[:, 0], 'blue')
axes[1].set_ylabel(r'$\theta(t)$ [deg]')

axes[2].plot(time, np.rad2deg(x_non_lin_traj[:, 1]), '.k', label='Measured')
axes[2].plot(time, np.rad2deg(sol[:-1].reshape(2, num_nodes).T)[:, 1], 'blue', label='Simulated')
axes[2].set_ylabel(r'$\omega(t)$ [deg/s]')
axes[2].set_xlabel('Time [s]')

leg = plt.legend()
```

Here I check if the constraints are satisfied:


```
constraint_violations = prob.con(sol).reshape(2, prob.num_constraints / 2).T
plt.plot(constraint_violations)
leg = plt.legend([sm.latex(s, mode='inline') for s in x])
```

## Linear Model

Now I use the linear identification model and the data from the non-linear simulation to see how well it does.


```
prob = Problem(obj,
               obj_grad,
               xdot - linear_sym_id_rhs,
               x,
               num_nodes,
               interval,
               known_parameter_map=known_par_map,
               known_trajectory_map={a: acceleration},
               integration_method='midpoint')
```


```
prob.addOption('linear_solver', 'ma57')
```


```
initial_guess = np.hstack((x_non_lin_traj.T.flatten(), par_map[k]))
```


```
sol, info = prob.solve(initial_guess)
sm.Eq(k_id, sol[-1])
```


```
fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8, 8))

axes[0].plot(time, acceleration)
axes[0].set_ylabel(r'a(t) [$ms^{-2}$]')

axes[1].plot(time, np.rad2deg(x_non_lin_traj[:, 0]), '.k')
axes[1].plot(time, np.rad2deg(sol[:-1].reshape(2, num_nodes).T)[:, 0], 'blue')
axes[1].set_ylabel(r'$\theta(t)$ [deg]')

axes[2].plot(time, np.rad2deg(x_non_lin_traj[:, 1]), '.k', label='Measured')
axes[2].plot(time, np.rad2deg(sol[:-1].reshape(2, num_nodes).T)[:, 1], 'blue', label='Simulated')
axes[2].set_ylabel(r'$\omega(t)$ [deg/s]')
axes[2].set_xlabel('Time [s]')

leg = plt.legend()
```


```
constraint_violations = prob.con(sol).reshape(2, prob.num_constraints / 2).T
plt.plot(constraint_violations)
```

# Conclusion

It seems that it is impossible to identify the exact gain value if any process noise is present, i.e. the identification model differs from the true system. Cranking up the number of nodes does not reduce the error in the identified parameter. The only way to get the exact gain is to set the process noise to zero, or at least make the ratio of $\frac{a}{w}$ large. It is also impossible for the simpler identification models to track the trajectory of the true system due to the unmodeled process noise. There is a systematic error and this bias potentially exists in both direct and indirect indentification methods of closed loop systems. In the case of identifying a human controller in a motor control task, we have little to no idea of what the process noise could be, so it is difficult to design a proper perturbation level.

3d726224e8f740548880c1bcac1d235afba0439c
ipynb
Jupyter Notebook
notebooks/systematic_error.ipynb
csu-hmc/inverted-pendulum-sys-id-paper
c41d9f7db16a05bc2e6194b62d303908ff45c2ec
CC-BY-4.0
2020-09-08T13:07:01.000Z
2020-09-08T13:07:01.000Z
notebooks/systematic_error.ipynb
csu-hmc/inverted-pendulum-sys-id-paper
c41d9f7db16a05bc2e6194b62d303908ff45c2ec
CC-BY-4.0
2015-02-05T18:11:28.000Z
2015-03-20T19:02:29.000Z
notebooks/systematic_error.ipynb
csu-hmc/inverted-pendulum-sys-id-paper
c41d9f7db16a05bc2e6194b62d303908ff45c2ec
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES