## Dynamic Estimation Tuning

## Main.EstimatorTuning History

Show minor edits - Show changes to markup

Attach:mhe_tuning.mp4 (:html:) <video width="450" controls autoplay loop>

<source src="/do/uploads/Main/mhe_tuning.mp4" type="video/mp4"> Your browser does not support the video tag.

</video> (:htmlend:)

n = 1 # process model order

p = GEKKO()

- Change to True for MacOS

rmt = False p = GEKKO(remote=rmt)

- p.u = p.MV()

p.u = p.MV(value=2)

p.u = p.MV(value=0)

- variable
- p.y = p.SV() #measurement

p.y = p.SV(value=0) #measurement

- Intermediate

p.x = [p.Intermediate(p.u)]

- Variables

p.x.extend([p.Var() for _ in range(n)]) #state variables p.y = p.SV() #measurement

p.Equation(p.tau * p.y.dt() == -p.y + p.K * p.u)

p.Equations([p.tau/n * p.x[i+1].dt() == -p.x[i+1] + p.x[i] for i in range(n)]) p.Equation(p.y == p.K * p.x[n])

- MHE Model

m = GEKKO()

m.time = np.linspace(0,20,41) #0-20 by 0.5 -- discretization must match simulation

- Model

m = GEKKO(remote=rmt)

- 0-20 by 0.5 -- discretization must match simulation

m.time = np.linspace(0,20,41)

- m.K = m.FV(value=3, lb=1, ub=3) #gain

m.K = m.FV(value=1, lb=1, ub=3) #gain

- m.tau = m.FV(value=4, lb=1, ub=10) #time constant

m.K = m.FV(value=1, lb=0.3, ub=3) #gain

- m.y = m.CV() #measurement

m.y = m.CV(value=0) #measurement

m.x = m.SV() #state variable m.y = m.CV() #measurement

m.Equation(m.tau * m.y.dt() == -m.y + m.K*m.u)

m.Equations([m.tau * m.x.dt() == -m.x + m.u,

m.y == m.K * m.x])

m.options.DIAGLEVEL = 0

m.y.STATUS = 1

m.K.DMAX = 0.0

- m.tau.DMAX = .1

m.tau.DMAX = 0.0

m.K.DMAX = 1 m.tau.DMAX = .1

m.y.MEAS_GAP = 0.25

m.y.TR_INIT = 1

m.y.MEAS_GAP = 0.0

- time vector

tm = np.linspace(0,25,51)

- noise = 0.25

noise = 0.0

- run process, estimator and control for cycles

y_meas = np.empty(cycles)

noise = 0.25

- values of u change randomly over time every 10th step

u_meas = np.zeros(cycles) step_u = 0 for i in range(0,cycles):

if ((i-1)%10) == 0: # random step (-5 to 5) step_u = step_u + (random()-0.5)*10 u_meas[i] = step_u

- run process and estimator for cycles

y_meas = np.zeros(cycles)

u_cont = np.empty(cycles) u = 2.0

- Create plot

plt.figure(figsize=(10,7)) plt.ion() plt.show()

# change input (u) if i==10: u = 3.0 elif i==20: u = 4.0 elif i==30: u = 1.0 elif i==40: u = 3.0 u_cont[i] = u ## process simulator #load u value p.u.MEAS = u_cont[i] #simulate

# process simulator p.u.MEAS = u_meas[i]

#load output with white noise y_meas[i] = p.y.MODEL + (random()-0.5)*noise ## estimator #load input and measured output m.u.MEAS = u_cont[i]

r = (random()-0.5)*noise y_meas[i] = p.y.value[1] + r # add noise # estimator m.u.MEAS = u_meas[i]

#optimize parameters

#store results y_est[i+1] = m.y.MODEL k_est[i+1] = m.K.NEWVAL tau_est[i+1] = m.tau.NEWVAL

y_est[i] = m.y.MODEL k_est[i] = m.K.NEWVAL tau_est[i] = m.tau.NEWVAL

plt.plot(y_meas[0:i],'b-') plt.plot(y_est[0:i],'r--')

plt.plot(tm[0:i+1],y_meas[0:i+1],'b-') plt.plot(tm[0:i+1],y_est[0:i+1],'r--')

plt.plot(np.ones(i)*p.K.value[0],'b-') plt.plot(k_est[0:i],'r--')

plt.plot(tm[0:i+1],np.ones(i+1)*p.K.value[0],'b-') plt.plot(tm[0:i+1],k_est[0:i+1],'r--')

plt.plot(np.ones(i)*p.tau.value[0],'b-') plt.plot(tau_est[0:i],'r--')

plt.plot(tm[0:i+1],np.ones(i+1)*p.tau.value[0],'b-') plt.plot(tm[0:i+1],tau_est[0:i+1],'r--')

plt.plot(u_cont[0:i])

plt.plot(tm[0:i+1],u_meas[0:i+1],'b-')

- Import packages

- Solve local or remote

rmt = True

n = 1 #process model order

p.u = p.MV()

- p.u = p.MV()

p.u = p.MV(value=2)

- Intermediate

p.x = [p.Intermediate(p.u)]

- Variables

p.x.extend([p.Var() for _ in range(n)]) #state variables p.y = p.SV() #measurement

- variable
- p.y = p.SV() #measurement

p.y = p.SV(value=0) #measurement

p.Equations([p.tau/n * p.x[i+1].dt() == -p.x[i+1] + p.x[i] for i in range(n)]) p.Equation(p.y == p.K * p.x[n])

p.Equation(p.tau * p.y.dt() == -p.y + p.K * p.u)

- p.u.FSTATUS = 1
- p.u.STATUS = 0
- Model

- MHE Model

- m.K = m.FV(value=3, lb=1, ub=3) #gain

- m.tau = m.FV(value=4, lb=1, ub=10) #time constant

m.x = m.SV() #state variable m.y = m.CV() #measurement

- m.y = m.CV() #measurement

m.y = m.CV(value=0) #measurement

m.Equations([m.tau * m.x.dt() == -m.x + m.u,

m.y == m.K * m.x])

m.Equation(m.tau * m.y.dt() == -m.y + m.K*m.u)

m.options.DIAGLEVEL = 0

m.K.DMAX = 1 m.tau.DMAX = .1

m.K.DMAX = 0.0

- m.tau.DMAX = .1

m.tau.DMAX = 0.0

noise = 0.25

- values of u change randomly over time every 10th step

u_meas = np.zeros(cycles) step_u = 0 for i in range(0,cycles):

if (i % 10) == 0: # random step (-5 to 5) step_u = step_u + (random()-0.5)*10 u_meas[i] = step_u

- run process and estimator for cycles

- noise = 0.25

noise = 0.0

- run process, estimator and control for cycles

y_est = np.empty(cycles) k_est = np.empty(cycles) tau_est = np.empty(cycles) for i in range(cycles):

# process simulator p.u.MEAS = u_meas[i] p.solve(remote=rmt)

y_est = np.zeros(cycles) k_est = np.zeros(cycles)*np.nan tau_est = np.zeros(cycles)*np.nan u_cont = np.empty(cycles) u = 2.0

- Create plot

plt.figure(figsize=(10,7)) plt.ion() plt.show()

for i in range(cycles-1):

# change input (u) if i==10: u = 3.0 elif i==20: u = 4.0 elif i==30: u = 1.0 elif i==40: u = 3.0 u_cont[i] = u ## process simulator #load u value p.u.MEAS = u_cont[i] #simulate p.solve() #load output with white noise

# estimator m.u.MEAS = u_meas[i]

## estimator #load input and measured output m.u.MEAS = u_cont[i]

m.solve(remote=rmt) y_est[i] = m.y.MODEL k_est[i] = m.K.NEWVAL tau_est[i] = m.tau.NEWVAL

- Plot results

plt.figure() plt.subplot(4,1,1) plt.plot(y_meas) plt.plot(y_est) plt.legend(('meas','pred')) plt.subplot(4,1,2) plt.plot(np.ones(cycles)*p.K.value[0]) plt.plot(k_est) plt.legend(('actual','pred')) plt.subplot(4,1,3) plt.plot(np.ones(cycles)*p.tau.value[0]) plt.plot(tau_est) plt.legend(('actual','pred')) plt.subplot(4,1,4) plt.plot(u_meas) plt.legend('u') plt.show()

#optimize parameters m.solve() #store results y_est[i+1] = m.y.MODEL k_est[i+1] = m.K.NEWVAL tau_est[i+1] = m.tau.NEWVAL plt.clf() plt.subplot(4,1,1) plt.plot(y_meas[0:i],'b-') plt.plot(y_est[0:i],'r--') plt.legend(('meas','pred')) plt.ylabel('y') plt.subplot(4,1,2) plt.plot(np.ones(i)*p.K.value[0],'b-') plt.plot(k_est[0:i],'r--') plt.legend(('actual','pred')) plt.ylabel('k') plt.subplot(4,1,3) plt.plot(np.ones(i)*p.tau.value[0],'b-') plt.plot(tau_est[0:i],'r--') plt.legend(('actual','pred')) plt.ylabel('tau') plt.subplot(4,1,4) plt.plot(u_cont[0:i]) plt.legend('u') plt.draw() plt.pause(0.05)

Tuning typically involves adjustment of objective function terms or constraints that limit the rate of change (DMAX), penalize the rate of change (DCOST), or set absolute bounds (LOWER and UPPER). Measurement availability is indicated by the parameter (FSTATUS). The optimizer can also include (1=on) or exclude (0=off) a certain adjustable parameter (FV) or manipulated variable (MV) with STATUS. Another important tuning consideration is the time horizon length. Including more points in the time horizon allows the estimator to reconcile the model to more data but also increases computational time. Below are common application, FV, MV, and CV tuning constants that are adjusted to achieve desired model predictive control performance.

Tuning typically involves adjustment of objective function terms or constraints that limit the rate of change (DMAX), penalize the rate of change (DCOST), or set absolute bounds (LOWER and UPPER). Measurement availability is indicated by the parameter (FSTATUS). The optimizer can also include (1=on) or exclude (0=off) a certain adjustable parameter (FV) or manipulated variable (MV) with STATUS. Another important tuning consideration is the time horizon length. Including more points in the time horizon allows the estimator to reconcile the model to more data but also increases computational time.

Below are common application, FV, MV, and CV tuning constants that are adjusted to achieve desired model predictive control performance.

Application constants are modified by indicating that the constant belongs to the group *nlc*. IMODE is adjusted to either solve the MHE problem with a simultaneous (5) or sequential (8) method. In the case below, the application IMODE is changed to sequential mode.

apm_option(server,app,'nlc.IMODE',8)

Application constants are modified by indicating that the constant belongs to the group *apm*. IMODE is adjusted to either solve the MHE problem with a simultaneous (5) or sequential (8) method. In the case below, the application IMODE is changed to simultaneous mode.

apm_option(server,app,'apm.IMODE',5)

**Objective:** Design an estimator to predict an unknown parameters so that a simple model is able to predict the response of a more complex process. Tune the estimator to achieve either tracking or predictive performance. *Estimated time: 2 hours.*

**Objective:** Design an estimator to predict an unknown parameters so that a simple model is able to predict the response of a more complex process. Tune the estimator to achieve either tracking or predictive performance. *Estimated time: 1 hour.*

- Manipulated Variable (MV) tuning

- Fixed Value (FV) - single parameter value over time horizon
- DMAX = maximum that FV can move each cycle
- LOWER = lower FV bound
- FSTATUS = feedback status with 1=measured, 0=off
- STATUS = turn on (1) or off (0) FV
- UPPER = upper FV bound

- Manipulated Variable (MV) - parameter can change over time horizon

(:toggle hide gekko button show="Show GEKKO (Python) Code":) (:div id=gekko:) (:source lang=python:)

- Import packages

import numpy as np from random import random from gekko import GEKKO import matplotlib.pyplot as plt

- Solve local or remote

rmt = True

- Process

p = GEKKO()

p.time = [0,.5]

n = 1 #process model order

- Parameters

p.u = p.MV() p.K = p.Param(value=1) #gain p.tau = p.Param(value=5) #time constant

- Intermediate

p.x = [p.Intermediate(p.u)]

- Variables

p.x.extend([p.Var() for _ in range(n)]) #state variables p.y = p.SV() #measurement

- Equations

p.Equations([p.tau/n * p.x[i+1].dt() == -p.x[i+1] + p.x[i] for i in range(n)]) p.Equation(p.y == p.K * p.x[n])

- options

p.options.IMODE = 4

- p.u.FSTATUS = 1
- p.u.STATUS = 0
- Model

m = GEKKO()

m.time = np.linspace(0,20,41) #0-20 by 0.5 -- discretization must match simulation

- Parameters

m.u = m.MV() #input m.K = m.FV(value=1, lb=1, ub=3) #gain m.tau = m.FV(value=5, lb=1, ub=10) #time constant

- Variables

m.x = m.SV() #state variable m.y = m.CV() #measurement

- Equations

m.Equations([m.tau * m.x.dt() == -m.x + m.u,

m.y == m.K * m.x])

- Options

m.options.IMODE = 5 #MHE m.options.EV_TYPE = 1

- STATUS = 0, optimizer doesn't adjust value
- STATUS = 1, optimizer can adjust

m.u.STATUS = 0 m.K.STATUS = 1 m.tau.STATUS = 1 m.y.STATUS = 1

- FSTATUS = 0, no measurement
- FSTATUS = 1, measurement used to update model

m.u.FSTATUS = 1 m.K.FSTATUS = 0 m.tau.FSTATUS = 0 m.y.FSTATUS = 1

- DMAX = maximum movement each cycle

m.K.DMAX = 1 m.tau.DMAX = .1

- MEAS_GAP = dead-band for measurement / model mismatch

m.y.MEAS_GAP = 0.25

m.y.TR_INIT = 1

- problem configuration
- number of cycles

cycles = 50

- noise level

noise = 0.25

- values of u change randomly over time every 10th step

u_meas = np.zeros(cycles) step_u = 0 for i in range(0,cycles):

if (i % 10) == 0: # random step (-5 to 5) step_u = step_u + (random()-0.5)*10 u_meas[i] = step_u

- run process and estimator for cycles

y_meas = np.empty(cycles) y_est = np.empty(cycles) k_est = np.empty(cycles) tau_est = np.empty(cycles) for i in range(cycles):

# process simulator p.u.MEAS = u_meas[i] p.solve(remote=rmt) y_meas[i] = p.y.MODEL + (random()-0.5)*noise # estimator m.u.MEAS = u_meas[i] m.y.MEAS = y_meas[i] m.solve(remote=rmt) y_est[i] = m.y.MODEL k_est[i] = m.K.NEWVAL tau_est[i] = m.tau.NEWVAL

- Plot results

plt.figure() plt.subplot(4,1,1) plt.plot(y_meas) plt.plot(y_est) plt.legend(('meas','pred')) plt.subplot(4,1,2) plt.plot(np.ones(cycles)*p.K.value[0]) plt.plot(k_est) plt.legend(('actual','pred')) plt.subplot(4,1,3) plt.plot(np.ones(cycles)*p.tau.value[0]) plt.plot(tau_est) plt.legend(('actual','pred')) plt.subplot(4,1,4) plt.plot(u_meas) plt.legend('u') plt.show() (:sourceend:) (:divend:)

- SOLVER
- 0=Try all available solvers
- 1=APOPT (MINLP, Active Set SQP)
- 2=BPOPT (NLP, Interior Point, SQP)
- 3=IPOPT (NLP, Interior Point, SQP)

- SOLVER

apm_option(server,app,'y.MEAS_GAS',1.0)

apm_option(server,app,'y.MEAS_GAP',1.0)

<iframe width="560" height="315" src="https://www.youtube.com/embed/yQWgSByYjd8?rel=0" frameborder="0" allowfullscreen></iframe>

<iframe width="560" height="315" src="https://www.youtube.com/embed/yw_a9ektOqc?rel=0" frameborder="0" allowfullscreen></iframe>

- CV_TYPE = CV type with 1=l_1-norm, 2=squared error

- SP = set point with CV_TYPE = 2
- SPLO = lower set point with CV_TYPE = 1
- SPHI = upper set point with CV_TYPE = 1
- STATUS = turn on (1) or off (0) MV
- TAU = reference trajectory time-constant
- TR_INIT = trajectory type, 0=dead-band, 1,2=trajectory
- TR_OPEN = opening at initial point of trajectory compared to end

- MEAS_GAP = measurement gap for estimator dead-band

Tuning typically involves adjustment of objective function terms or constraints that limit the rate of change (DMAX), penalize the rate of change (DCOST), or set absolute bounds (LOWER and UPPER). Measurement availability is indicated by the parameter (FSTATUS). The optimizer can also include (1=on) or exclude (0=off) a certain adjustable parameter (FV) or manipulated variable (MV) with STATUS. Below are common application, FV, MV, and CV tuning constants that are adjusted to achieve desired model predictive control performance.

Tuning typically involves adjustment of objective function terms or constraints that limit the rate of change (DMAX), penalize the rate of change (DCOST), or set absolute bounds (LOWER and UPPER). Measurement availability is indicated by the parameter (FSTATUS). The optimizer can also include (1=on) or exclude (0=off) a certain adjustable parameter (FV) or manipulated variable (MV) with STATUS. Another important tuning consideration is the time horizon length. Including more points in the time horizon allows the estimator to reconcile the model to more data but also increases computational time. Below are common application, FV, MV, and CV tuning constants that are adjusted to achieve desired model predictive control performance.

- IMODE = 5 or 8 for moving horizon estimation

- IMODE = 5 or 8 for moving horizon estimation

apm_option(server,app,'nlc.IMODE',8)

apm_option(server,app,'nlc.IMODE',8)

#### Exercise

**Objective:** Design an estimator to predict an unknown parameters so that a simple model is able to predict the response of a more complex process. Tune the estimator to achieve either tracking or predictive performance. *Estimated time: 2 hours.*

Design an estimator to predict *K* and *tau* of a 1st order model to predict the dynamic response of a 1st order, 2nd order, and 10th order process. For the 2nd and 10th order systems, there is process/model mismatch. This means that the structure of the model can never exactly match the actual process because the equations are inherently incorrect. The parameter values are adjusted to best approximate the process even though the model is deficient. The process order is adjusted in the file *process.apm* file in the *Constants* section.

Constants ! process model order n = 1 ! change to 1, 2, and 10

In each case, tune the estimator to favor either acceptable tracking or predictive performance. Tracking performance is the ability of the estimator to synchronize with measurements and is demonstrated with overall agreement between the model predictions and the measurements. Predictive performance sacrifices tracking performance to achieve more consistent values that are valid over a longer predictive horizon for model predictive control.

#### Solution

(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/yQWgSByYjd8?rel=0" frameborder="0" allowfullscreen></iframe> (:htmlend:)

The dynamic estimation tuning is the process of adjusting certain objective function terms to give more desirable solutions. As an example, dynamic estimation application may either track noisy data too closely or the updates may be too slow to catch unmeasured disturbances of interest.

Dynamic estimation tuning is the process of adjusting certain objective function terms to give more desirable solutions. As an example, a dynamic estimation application such as moving horizon estimation (MHE) may either track noisy data too closely or the updates may be too slow to catch unmeasured disturbances of interest. Tuning is the process of achieving acceptable estimator performance based on unique aspects of the application.

#### Common Tuning Parameters for MHE

Tuning typically involves adjustment of objective function terms or constraints that limit the rate of change (DMAX), penalize the rate of change (DCOST), or set absolute bounds (LOWER and UPPER). Measurement availability is indicated by the parameter (FSTATUS). The optimizer can also include (1=on) or exclude (0=off) a certain adjustable parameter (FV) or manipulated variable (MV) with STATUS. Below are common application, FV, MV, and CV tuning constants that are adjusted to achieve desired model predictive control performance.

- Application tuning
- IMODE = 5 or 8 for moving horizon estimation
- DIAGLEVEL = diagnostic level (0-10) for solution information
- EV_TYPE = 1 for l
_{1}-norm and 2 for squared error objective - MAX_ITER = maximum iterations
- MAX_TIME = maximum time before stopping
- MV_TYPE = Set default MV type with 0=zero-order hold, 1=linear interpolation

- Manipulated Variable (MV) tuning
- COST = (+) minimize MV, (-) maximize MV
- DCOST = penalty for MV movement
- DMAX = maximum that MV can move each cycle
- FSTATUS = feedback status with 1=measured, 0=off
- LOWER = lower MV bound
- MV_TYPE = MV type with 0=zero-order hold, 1=linear interpolation
- STATUS = turn on (1) or off (0) MV
- UPPER = upper MV bound

- Controlled Variable (CV) tuning
- COST = (+) minimize MV, (-) maximize MV
- CV_TYPE = CV type with 1=l_1-norm, 2=squared error
- FSTATUS = feedback status with 1=measured, 0=off
- SP = set point with CV_TYPE = 2
- SPLO = lower set point with CV_TYPE = 1
- SPHI = upper set point with CV_TYPE = 1
- STATUS = turn on (1) or off (0) MV
- TAU = reference trajectory time-constant
- TR_INIT = trajectory type, 0=dead-band, 1,2=trajectory
- TR_OPEN = opening at initial point of trajectory compared to end

There are several ways to change the tuning values. Tuning values can either be specified before an application is initialized or while an application is running. To change a tuning value before the application is loaded, use the *apm_option()* function such as the following example to change the lower bound in MATLAB or Python for the FV named *p*.

apm_option(server,app,'p.LOWER',0)

The upper and lower measurement deadband for a CV named *y* are set to values around the measurement. In this case, an acceptable range for the model prediction is to intersect the measurement of 10.0 between 9.5 and 10.5 with a MEAS_GAP of 1.0 (or +/-0.5).

apm_option(server,app,'y.MEAS_GAS',1.0)

Application constants are modified by indicating that the constant belongs to the group *nlc*. IMODE is adjusted to either solve the MHE problem with a simultaneous (5) or sequential (8) method. In the case below, the application IMODE is changed to sequential mode.

apm_option(server,app,'nlc.IMODE',8)

(:title Dynamic Estimation Tuning:) (:keywords tuning, Python, MATLAB, Simulink, moving horizon, time window, dynamic data, validation, estimation, differential, algebraic, tutorial:) (:description Tuning of an estimator for improved rejection of corrupted data with outliers, drift, and noise:)

The dynamic estimation tuning is the process of adjusting certain objective function terms to give more desirable solutions. As an example, dynamic estimation application may either track noisy data too closely or the updates may be too slow to catch unmeasured disturbances of interest.