ESRF Double Crystal Monochromator - Dynamical Multi-Body Model

Let’s considered the system \(\bm{G}(s)\) with:

• 3 inputs: force applied to the 3 fast jacks
• 3 outputs: measured displacement by the 3 interferometers pointing at the ring second crystal

It is schematically shown in Figure 1.

The system is identified from the Simscape model.

%% Input/Output definition
clear io; io_i = 1;

%% Inputs
% Control Input {3x1} [N]
io(io_i) = linio([mdl, '/control_system'], 1, 'openinput');  io_i = io_i + 1;

%% Outputs
% Interferometers {3x1} [m]
io(io_i) = linio([mdl, '/DCM'], 1, 'openoutput'); io_i = io_i + 1;

%% Extraction of the dynamics
G = linearize(mdl, io);

size(G)

size(G)
State-space model with 3 outputs, 3 inputs, and 24 states.

load('dcm_kinematics.mat');

Using the forward and inverse kinematics, we can computed the dynamics from piezo forces to axial motion of the 3 fastjacks (see Figure 2).

%% Compute the system in the frame of the fastjacks
G_pz = J_a_h*inv(J_2h_s)*G;

The DC gain of the new system shows that the system is well decoupled at low frequency.

dcgain(G_pz)

The bode plot of \(\bm{G}_{\text{fj}}(s)\) is shown in Figure 3.

G_pz = diag(1./diag(dcgain(G_pz)))*G_pz;

G_mr = inv(J_s_r)*G*inv(J_a_r');

dcgain(G_mr)

This results in a coupled system. The main reason is that, as we map forces to the center of the ring crystal and not at the center of mass/stiffness of the moving part, vertical forces will induce rotation and torques will induce vertical motion.

Interferometer noise:

Wn = 6e-11*(1 + s/2/pi/200)/(1 + s/2/pi/60); % m/sqrt(Hz)

Measurement noise: 0.79 [nm,rms]

DAC noise (amplified by the PI voltage amplifier, and converted to newtons):

Wdac = tf(3e-8); % V/sqrt(Hz)
Wu = Wdac*22.5*10; % N/sqrt(Hz)

DAC noise: 0.95 [uV,rms]

Disturbances:

Wd = 5e-7/(1 + s/2/pi); % m/sqrt(Hz)

Disturbance motion: 0.61 [um,rms]

%% Save ASD of noise and disturbances
save('mat/asd_noises_disturbances.mat', 'Wn', 'Wu', 'Wd');

The comparison of the amplitude spectral density of the measurement noise and of the jack parasitic motion is performed in Figure 6. It confirms that the sensor noise is low enough to measure the motion errors of the crystal.

Krad_g1 = eye(3)*s/(s^2/(2*pi*500)^2 + 2*s/(2*pi*500) + 1);

As can be seen in Figure 8, very little damping can be added using relative damping strategy using strain gauges.

Krad = -g*Krad_g1;

The controller is implemented on Simscape, and the damped plant is identified.

%% Input/Output definition
clear io; io_i = 1;

%% Inputs
% Control Input {3x1} [N]
io(io_i) = linio([mdl, '/control_system'], 1, 'input');  io_i = io_i + 1;

%% Outputs
% Force Sensor {3x1} [m]
io(io_i) = linio([mdl, '/DCM'], 1, 'openoutput'); io_i = io_i + 1;

%% DCM Kinematics
load('dcm_kinematics.mat');

%% Identification of the Open Loop plant
controller.type = 0; % Open Loop
G_ol = J_a_r*inv(J_s_r)*linearize(mdl, io);
G_ol.InputName  = {'u_ur',  'u_uh',  'u_d'};
G_ol.OutputName = {'d_ur',  'd_uh',  'd_d'};

%% Identification of the damped plant with Relative Active Damping
controller.type = 2; % RAD
G_dp = J_a_r*inv(J_s_r)*linearize(mdl, io);
G_dp.InputName  = {'u_ur',  'u_uh',  'u_d'};
G_dp.OutputName = {'d_ur',  'd_uh',  'd_d'};

%% Input/Output definition
clear io; io_i = 1;

%% Inputs
% Control Input {3x1} [N]
io(io_i) = linio([mdl, '/control_system'], 1, 'openinput');  io_i = io_i + 1;

%% Outputs
% Force Sensor {3x1} [m]
io(io_i) = linio([mdl, '/DCM'], 3, 'openoutput'); io_i = io_i + 1;

%% Extraction of the dynamics
G_fs = linearize(mdl, io);

The Bode plot of the identified dynamics is shown in Figure 10. At high frequency, the diagonal terms are constants while the off-diagonal terms have some roll-off.

We want to have integral action around the resonances of the system, but we do not want to integrate at low frequency. Therefore, we can use a low pass filter.

%% Integral Force Feedback Controller
Kiff_g1 = eye(3)*1/(1 + s/2/pi/20);

%% Integral Force Feedback Controller with optimal gain
Kiff = g*Kiff_g1;

%% Save the IFF controller
save('mat/Kiff.mat', 'Kiff');

Both the Open Loop dynamics (see Figure 2) and the dynamics with IFF (see Figure 12) are identified.

We are here interested in the dynamics from \(\bm{u}^\prime = [u_{u_r}^\prime,\ u_{u_h}^\prime,\ u_d^\prime]\) (input of the damped plant) to \(\bm{d}_{\text{fj}} = [d_{u_r},\ d_{u_h},\ d_d]\) (motion of the crystal expressed in the frame of the fast-jacks). This is schematically represented in Figure 12.

The dynamics from \(\bm{u}\) to \(\bm{d}_{\text{fj}}\) (open-loop dynamics) and from \(\bm{u}^\prime\) to \(\bm{d}_{\text{fs}}\) are compared in Figure 13. It is clear that the Integral Force Feedback control strategy is very effective in damping the resonances of the plant.

Let’s identify the damped plant.

Let’s design a controller with a bandwidth of 100Hz. As the plant is well decoupled and well approximated by a constant at low frequency, the high authority controller can easily be designed with SISO loop shaping.

%% Controller design
wc = 2*pi*100; % Wanted crossover frequency [rad/s]
a = 2; % Lead parameter

Khac = diag(1./diag(abs(evalfr(G_dp, 1j*wc)))) * ... % Diagonal controller
        wc/s * ... % Integrator
        1/(sqrt(a))*(1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a))) * ... % Lead
        1/(s^2/(4*wc)^2 + 2*s/(4*wc) + 1); % Low pass filter

%% Save the HAC controller
save('mat/Khac_iff.mat', 'Khac');

%% Loop Gain
L_hac_lac = G_dp * Khac;

As shown in the Root Locus plot in Figure 17, the closed loop system should be stable.

In order to estimate the performances of the HAC-IFF control strategy, the transfer function from motion errors of the stepper motors to the motion error of the crystal is identified both in open loop and with the HAC-IFF strategy.

It is first verified that the closed-loop system is stable:

isstable(T_hl)

1

And both transmissibilities are compared in Figure 18.

Let’s compute the amplitude spectral density of the jack motion errors due to the sensor noise, the actuator noise and disturbances.

%% Computation of ASD of contribution of inputs to the closed-loop motion
% Error due to disturbances
asd_d = abs(squeeze(freqresp(Wd*(1/(1 + G_dp(1,1)*Khac(1,1))), f, 'Hz')));
% Error due to actuator noise
asd_u = abs(squeeze(freqresp(Wu*(G_dp(1,1)/(1 + G_dp(1,1)*Khac(1,1))), f, 'Hz')));
% Error due to sensor noise
asd_n = abs(squeeze(freqresp(Wn*(G_dp(1,1)*Khac(1,1)/(1 + G_dp(1,1)*Khac(1,1))), f, 'Hz')));

The closed-loop ASD is then:

%% ASD of the closed-loop motion
asd_cl = sqrt(asd_d.^2 + asd_u.^2 + asd_n.^2);

The obtained ASD are shown in Figure 19.

Let’s compare the open-loop and close-loop cases (Figure 20).