Overview

This project was completed individually for 4CM120 — Data-Based Optimization of Control Systems (TU Eindhoven, MSc Mechanical Engineering, Feb–Apr 2026). Taught by Dr.ir. Thijs van Keulen and Dr.ir. Leroy Hazeleger, the course studies real-time, model-free optimization for high-tech mechatronic systems.

The problem — nanometer-scale localization of a moving stage with respect to a reference frame — is directly relevant to atomic force microscopes and lithography machines, where fast and accurate alignment is critical. The approach implemented here is based on the patent by Van Keulen et al. (2021, Method of controlling a position of a first object relative to a second object, PCT/EP2020/070979).

The objective function is the light intensity measured by a photodiode on the moving body as it passes under a fixed reference lens:

This Gaussian has a unique maximum at $x^* = 120\,\text{rad}$ with $J^* = 0.554$. The goal of the ESC controller is to find this peak without any analytical model of $Q_J$ — using only the measured intensity signal $J(t)$.

The Assignment

The project is split into two deliverables — both embedded below:

The Plant Model

The moving stage uses a combination of feedforward and feedback control to follow a reference position $r$. The underlying fourth-order transfer function from input force $u$ to position $x$ is:

with a PD + second-order Low-Pass feedback controller:

where $k_p = 0.05$, $k_d = 0.005$, $\omega_\text{lp} = 400\pi\,\text{rad/s}$, $\zeta_\text{lp} = 0.7$. The full loop runs discrete-time at $F_s = 4\,\text{kHz}$ via zero-order-hold discretization.

Part I — Classical ESC in Simulation (Assignment 1)

How Extremum Seeking Control Works

ESC estimates the gradient of the unknown objective $Q_J$ by injecting a sinusoidal dither signal $a\cos(\omega t)$ onto the reference. The measured output is then demodulated (multiplied by the same cosine) and filtered to extract the gradient component $\xi(t)$. An integrator uses this estimate to steer the nominal position $\hat{x}$ toward the optimum.

The loop has five elements:

1. Dither signal — $a\cos(\omega t)$ added to $\hat{x}$
2. High-pass filter $H_\text{HP}(s) = s/(s + \omega_f)$ — removes DC
3. Demodulation — multiplication by $\frac{2}{a}\cos(\omega t)$ extracts gradient information
4. Low-pass filter $H_\text{LP}(s) = \omega_f/(s + \omega_f)$ — smooths the estimate
5. Integrator with gain $k$ — drives $\hat{x}$ toward increasing $J$

A critical design constraint is the timescale separation: $\omega_f \ll \omega \ll \omega_\text{BW}$.

Unimodality Assumption (Q1)

For ESC to converge to the global optimum from any initial condition $x_0$, the objective function must be unimodal — i.e., have a single critical point. Since $Q_J$ is Gaussian, $Q_J''(120) < 0$ and $|x - 120|$ causes strict monotonic decrease in both directions. Only one zero of the gradient exists at $x^* = 120\,\text{rad}$, so the assumption is satisfied.

Optimal Parameters on the Static Map (Q2)

Starting from $x_0 = 10\,\text{rad}$ with $\omega = 1\,\text{Hz}$ prescribed, tuning converged on:

• Convergence time: 17.3 s
• Steady-state error: < 0.1 rad
• Peak intensity reached: $J = 0.554$

Measured vs. Analytical Objective (Q3)

The analytical $Q_J$ was replaced with an experimentally measured objective function obtained by sweeping the stage slowly past the photodiode. Intensity values were binned into 0.5 rad intervals and averaged, producing a lookup table for Simulink interpolation.

Key differences:

• Peak offset: the experimental peak appears at $x \approx 98\,\text{rad}$ rather than $120\,\text{rad}$ — an encoder-zero offset
• Steeper slopes: the real lens-photodiode geometry produces sharper gradients near the peak
• The analytical parameters caused divergence on the measured objective because the high filter cutoff created phase lag through the steeper gradient

Re-tuning for the measured case required a lower filter cutoff ($\omega_f = 1.01\,\text{rad/s}$) and slightly adjusted gains ($k = 5250$, $a = 7\,\text{rad}$).

Lookup Table vs Analytical Function (Q3)

Incorporating Plant Dynamics (Q5–Q7)

So far the plant was treated as a static map. In reality the moving stage has non-trivial dynamics: for ESC to work, the dynamic system $G: r \to J$ must behave as a static map at the dither frequency — i.e., the closed-loop transfer function $T(s) = x(s)/r(s)$ must satisfy:

The Bode plot confirms a closed-loop bandwidth of $\omega_\text{BW} = 67.3\,\text{rad/s}$. At the 1 Hz dither frequency ($\omega = 6.28\,\text{rad/s}$), $|T(j\omega)| = -0.66\,\text{dB}$ and $\angle T(j\omega) = -9.3°$ — well within the static-map regime.

With dynamics and noise enabled (disturbances = true, dynamics = true), the ESC was re-calibrated:

Despite the new constraints (phase lag, stage resonances, tighter filter speeds), convergence was achieved in ~20.8 s with all requirements met.

Part II — ESC on Physical Hardware (Assignment 2)

Assignment 2 took the simulation framework and deployed it on the real test setup — a motorized rotary stage with a photodiode. Three distinct ESC architectures were compared, both in simulation and on hardware:

1. Classical ESC — sinusoidal dither with HP/LP filters
2. Moving Average ESC — replaces the HP/LP filters with an exact moving average
3. Sampled-Data ESC — discrete step perturbations with finite-difference gradients

Classical ESC on Hardware (C-Q1)

The Q7 simulation parameters were used as a starting point. Four experimental runs were performed on the physical setup.

The best run matched the Assignment 1 Q7 parameters exactly and converged in 9.1 s — over a second faster than the simulation (10.6 s), due to the steeper real gradient.

Encoder Offset Discovery

The physical encoder's zero position did not match the model's $x = 0$. The system converged to $\hat{r}_\text{ss} \approx 60\text{–}67\,\text{rad}$ rather than the model's 120 rad, but the peak intensities were identical ($J \approx 0.555$). Aligning the data with a ~53 rad offset confirmed the experimental shape follows the model Gaussian near the peak but with steeper slopes further out.

Simulation vs. Experiment

• Convergence speed: within 1.5 s of each other, validating the model
• Steady-state oscillations: smaller on the experiment ($\pm 0.03\,\text{rad}$ vs. $\pm 0.09\,\text{rad}$), suggesting additional unmodeled HF attenuation
• Gradient estimate $\xi(t)$: less volatile in simulation, hinting at unmodeled plant dynamics

Positioning-Accuracy Trade-off (C-Q2)

A subtle but important effect: the actual position is $r(t) = \hat{r}(t) + a\cos(\omega t)$, so even when $\hat{r}$ has converged exactly to $x^*$, the position still oscillates by $\pm a$ around the optimum. If $a > 0.1\,\text{rad}$, the 0.1 rad positioning requirement is not met just by the nominal position alone.

Three levers affect accuracy:

• Reduce $a$ — smaller dither reduces oscillation but weakens the gradient signal (noisier $\xi$, slower convergence)
• Reduce $k$ — less overshoot but slower convergence
• Reduce $\omega_f$ — better noise rejection but slower gradient extraction

This creates a fundamental convergence-speed vs. positioning-accuracy trade-off. For the dither-based accuracy requirement $|\hat{r} - x^*| + a < 0.1\,\text{rad}$, both $a < 0.1$ and $\hat{r}$ must converge tightly to $x^*$.

Moving Average ESC (C-Q3 & C-Q4)

The moving average (MA) filter replaces the HP → demodulation → LP chain with a single operation: averaging the demodulated signal over exactly one dither period $T = 2\pi/\omega$.

The theoretical benefits are significant:

• Exact gradient extraction — the MA over one cosine period is exactly zero, perfectly canceling the dither component and all harmonics
• No phase lag — the classical HP/LP filters introduce phase lag $-\arctan(\omega/\omega_f)$; the MA has only a fixed $T/2$ delay that is frequency-independent
• Simpler tuning — only the window length $T$ needs to match the dither period
• Higher dither frequencies possible — no strict $\omega_f \ll \omega$ timescale separation required

MA Simulation Results

By exploiting the relaxed timescale constraint, the dither frequency was increased from $\omega = 2\pi$ to $\omega = 24.5\,\text{rad/s}$. The MA window shortened from $T = 1.0\,\text{s}$ to $T = 0.26\,\text{s}$, enabling more aggressive tuning: $k = 24000$, $a = 15\,\text{rad}$. Simulation convergence: 3.9 s — a 5× improvement over the classical approach.

MA Experimental Results

Seven runs on the physical setup. The best run (MA2: $k = 35000$, $a = 5.8$, $f = 10\,\text{Hz}$) converged in 1.3 s — a 7× speedup over the classical experiment (9.1 s).

The MA approach proved sensitive to aggressive gain tuning: MA1 at $k = 45000$ went immediately unstable, showing the test setup has a lower maximum stable gain than simulation. Lower dither frequencies (2–5 Hz) converged significantly slower.

Classical vs. MA Head-to-Head

Sampled-Data ESC (SD-Q1 → SD-Q4)

The sampled-data approach replaces the continuous sinusoidal dither with discrete step perturbations of amplitude $\pm c_v$ and estimates the gradient via central differences:

where $J^\pm$ are steady-state intensity measurements at positions $\hat{r} \pm c_v$. A discrete-time integrator then updates: $\hat{r}_{k+1} = \hat{r}_k + \lambda \cdot \xi_k$.

Three parameters govern SD performance:

• Step amplitude $c_v$ — finite-difference accuracy vs. signal-to-noise ratio
• Waiting time $T$ — must exceed the plant's settling time for unbiased measurements
• Step size $\lambda$ — analogous to the integrator gain $k$ in classical ESC

Dither-Amplitude Analysis (SD-Q1)

• Small $c_v$ (< 0.5 rad): accurate finite difference, but $\Delta J$ is small relative to sensor noise → noisy estimates
• Medium $c_v$ (0.5–2 rad): good tracking, strong gradient signal
• Large $c_v$ (> 5 rad): overshoot + higher-order Taylor terms ($\propto c_v^2$) degrade the approximation

SD Simulation Results (SD-Q3)

Tuned parameters achieved convergence in 8.7 s:

SD Experimental Results (SD-Q4)

Twelve runs were performed on the physical setup. None converged within the 0.1 rad bound due to a discovered implementation bug: the discrete-time integrator gain was set to $\lambda$ instead of $\lambda/(2T)$. The best run (Run 8) still settled to within 3.1 rad of the optimum.

Three-Way Comparison

Key architectural differences:

• Classical: robust to noise (LP filter averages over many periods), but bound by $\omega_f \ll \omega \ll \omega_\text{BW}$
• Moving Average: exact dither cancellation enables higher frequencies and faster convergence, but sensitive to high gains and requires precise window matching
• Sampled-Data: simplest gradient estimation (only two measurements), naturally handles constraints via barrier functions, scales linearly to MIMO ($2N$ evaluations per iteration), but vulnerable to noise spikes on individual measurements

Engineering Takeaways

• Timescale separation ($\omega_f \ll \omega \ll \omega_\text{BW}$) is the fundamental design constraint for classical ESC — it limits the achievable convergence speed but ensures filter stability. The MA filter relaxes this constraint and delivers order-of-magnitude speedups.
• Dither amplitude creates a fundamental trade-off: larger $a$ means better gradient SNR but more steady-state oscillation around the optimum. For the 0.1 rad accuracy requirement, $a$ itself must be < 0.1 rad.
• Simulation-to-hardware transfer is non-trivial: encoder offsets, unmodeled dynamics, and lower stability margins all required parameter re-tuning between simulation and the physical test rig.
• Implementation details matter: the sampled-data integrator gain normalization $\lambda/(2T)$ vs. $\lambda$ was a subtle bug that made the difference between stable and unstable behavior on hardware.

Implementation Details

All three ESC approaches were deployed on the physical setup through MATLAB scripts that programmatically configure Simulink model parameters. The system sampling rate is $F_s = 4\,\text{kHz}$.

Classical / MA (ESC_CT.m)

The script toggles between classical and moving-average modes using gain switches, and discretizes the HP/LP filter prototypes before deployment:

% Classical ESC filter discretization
HPc = tf([1 0], [1 DE.highpass_filter_pole * 2 * pi]);
HPd = c2d(HPc, 1/Fs);
LPc = tf(1, [1 / (DE.highpass_filter_pole * 2 * pi) 1]);
LPd = c2d(LPc, 1/Fs);
 
% Toggle classical vs MA via gain switches
if DE.moving_average
    set_param(ModelName + "/DE/Gain_classic", "gain", "0")
    set_param(ModelName + "/DE/Gain_MA", "gain", "1")
end

For the MA filter, the FIR window length is $N = F_s \cdot 2\pi/\omega$ samples.

Sampled-Data (ESC_DT.m)

The discrete-time implementation uses a pulse generator with period $2T$ (one $+c_v$ phase and one $-c_v$ phase), triggered subsystems sample $J$ at the end of each phase, and a discrete-time integrator updates via forward Euler:

% Critical: integrator gain must account for sample time
set_param(ModelName + "/optimizer/Gain", "gain", ...
    num2str(optimizer.gain / (2 * dither.time)))

Objective-Function Measurement (linear_motion.m)

A linear sweep script ramps the stage position at constant velocity while recording encoder position and photodiode intensity. This produces a scatter plot of $(x, J)$ pairs that maps the lens-photodiode geometry — the source for the Q3 lookup table.

Technologies Used

MATLAB, Simulink, real-time hardware control, extremum seeking control, moving average filtering, sampled-data optimization, finite-difference gradient estimation, frequency-domain analysis (Bode plots, bandwidth), system identification, experimental parameter tuning