pole placement controller gain trouble with identified receptance

Hi,

I'm trying to determine the controller gains f and g that allow placing the desired poles p in a 2-DOF system. My code implements two methods:

• Method 1 (that gives the correct f and g) requires the receptance matrix H in "tf" data to do evalfr. However, when estimating H from real system measurements using tfest, the antiresonance near the resonance at around 6 Hz is not properly captured.
• Method 2 attempts to solve this issue by using tfestimate to obtain H_11 and H_21​ (in complex double instead of tf), which correctly show the antiresonance. Instead of evalfr, this method interpolates H_11​ and H_21​ around the target frequency f_target​, effectively replicating evalfr without needing a transfer function object because it only needs complex double data.

The issue is that Method 2 produces a gain f equal to zero, and the computed g value differs from Method 1. How can I modify the code to ensure that both methods yield equivalent gains f and g?

Thank you for your help!

P.S. The attached .mat file contains measurements from an idealized MATLAB model, which correctly shows the Bode plot with the expected antiresonance (then I will use the real model measures).

% Upload I/O data
load 'test_identification_analytical_model_.mat'

% Data measured from simulated model
th1abs = theta1_rad_MEAS.signals.values;   
th2abs = theta2_rad_MEAS.signals.values;
tauIn = tauInput.signals.values;

nfs = 4096*4;   % Number of samples for FFT
Fs = 1000;   % Sampling frequency

% Transfer functions estimate from input 1 (tauIn) to outputs 1 and 2 (th1abs and th2abs)
[H_11, f_H_11] = tfestimate(squeeze(tauIn), squeeze(th1abs), hann(floor(size(squeeze(tauIn), 1) / 10)), [], nfs*4, Fs);
[H_21, f_H_21] = tfestimate(squeeze(tauIn), squeeze(th2abs), hann(floor(size(squeeze(tauIn), 1) / 10)), [], nfs*4, Fs);

% Estimated tf smoothing, by movmean with [x,x] values
H_11 = movmean(H_11,[2,2]);
H_21 = movmean(H_21,[2,2]);

% Estimated tf plot
figure(1)
subplot(1,2,1);
semilogx(f_H_11, mag2db(abs(H_11)), 'b', 'LineWidth', 0.7);
hold on;
subplot(1,2,2);
semilogx(f_H_21, mag2db(abs(H_21)), 'b', 'LineWidth', 0.7);
hold on;

% Analytical tf plot
figure(2)
P = bodeoptions; 
P.FreqUnits = 'Hz';   % Frequency axys in Hz
bodeplot([H(1,1); H(2,1)], P);
grid on;
hold on;

%% POLE PLACEMENT

p = [-0.8925 + 1j*7.4483 ; -0.8925 - 1j*7.4483 ; -0.9204 + 1j*37.6548 ; -0.9204 - 1j*37.6548];   % System poles to place

B = [1;0];
n = size(B,1);
tau = 0;

%--------------------------------------------------
% 1st METHOD

G = zeros(2*n, 4); 

for i = 1:2*n
    r = evalfr(H, p(i))*B;   % to keep only H_11 and H_21
    G(i,:) = [r.' p(i)*r.'];
end

w = -G\(-exp(p*tau));

g = real(w(1:n))'

g = 1×2

1.0e-04 * -0.1471 0.0976

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

f = real(w(n+1:2*n))'

f = 1×2

1.0e-06 * -0.3688 -0.3196

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%--------------------------------------------------
% 2nd METHOD

G = zeros(2*n, 4); 

for i = 1:2*n
    f_target = abs(p(i)) / (2*pi);
    H1_interp = interp1(f_H_11, H_11, f_target, 'linear', 'extrap');
    H2_interp = interp1(f_H_21, H_21, f_target, 'linear', 'extrap');
    r = [H1_interp; H2_interp];  
    G(i,:) = [r.' p(i)*r.'];
end

w = -G\(-exp(p*tau));

g = real(w(1:n))'

g = 1×2

-0.0030 0.0052

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

f = real(w(n+1:2*n))'

f = 1×2

0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%--------------------------------------------------

[Sam Chak]

Hi @GIOVANNI

I am curious about your expectations of the said pole placement technique (which is commonly taught in textbooks, but generally not explained in detail). Based on the desired closed-loop poles, are you anticipating that the closed-loop system will exhibit such oscillatory behavior with a maximum overshoot of slightly under 69%?

%% desired closed-loop poles
p   = [-0.8925+1j*7.4483; -0.8925-1j*7.4483; -0.9204+1j*37.6548; -0.9204-1j*37.6548];   % System poles to place

%% idealized closed-loop system
G   = tf(zpk([], p, real(prod(p))))

G =
 
                     7.984e04
  ----------------------------------------------
  s^4 + 3.626 s^3 + 1478 s^2 + 2636 s + 7.984e04
 
Continuous-time transfer function.

step(G), grid on

stepinfo(G)

ans = struct with fields:
         RiseTime: 0.1321
    TransientTime: 4.2969
     SettlingTime: 4.2969
      SettlingMin: 0.5249
      SettlingMax: 1.6864
        Overshoot: 68.6374
       Undershoot: 0
             Peak: 1.6864
         PeakTime: 0.4188

[GIOVANNI]

The poles in my code are only for this example, because my problem is to fix the section that finds gains f and g.

If i place poles with higher real and imaginary part the dynamic will be faster and with less overshoot.

in addition to allocating the desired poles, I have a simulink scheme that uses an internal model control with a filter that calculates coefficients to introduce anti-resonances at the detected disturbance frequency.

However the system diverges if I do not use the correct gains f and g to place the poles and i can't do with evalfr with real model measures because there is a problem with tfest function that can't detect antiresonance in the estimated receptance H, so i'm trying to find f and g using output of tfestimate H_11 and H_21 (see figure 1 in the code)

[Sam Chak]

Hi @GIOVANNI

Edit: In my earlier post, I made some erroneous computations of the frequency in Hz instead of radians, which resulted in the incorrect display of the Bode diagrams. Here, I have made a few corrections to the frequency arguments in the frd() and bode() commands.

% Upload I/O data
load 'test_identification_analytical_model_.mat'
whos

  Name                 Size              Bytes  Class     Attributes

  H                    2x2                2350  tf                  
  ans                  1x54                108  char                
  tauInput             1x1             2400943  struct              
  theta1_rad_MEAS      1x1             2400943  struct              
  theta2_rad_MEAS      1x1             2400943  struct              

Htf     = tf(H.Numerator, H.Denominator);
H11tf   = Htf(1,1)

H11tf =
 
           113 s^2 + 201.6 s + 1.417e05
  ----------------------------------------------
  s^4 + 3.626 s^3 + 1478 s^2 + 2636 s + 7.984e04
 
Continuous-time transfer function.

H21tf   = Htf(2,1)

H21tf =
 
           -186 s^2 - 268.8 s + 1.31e05
  ----------------------------------------------
  s^4 + 3.626 s^3 + 1478 s^2 + 2636 s + 7.984e04
 
Continuous-time transfer function.

% Data measured from simulated model
th1abs = theta1_rad_MEAS.signals.values;   
th2abs = theta2_rad_MEAS.signals.values;
tauIn  = tauInput.signals.values;

nfs = 4096*4;   % Number of samples for FFT
Fs = 1000;   % Sampling frequency

% Transfer functions estimate from input 1 (tauIn) to outputs 1 and 2 (th1abs and th2abs)
[H_11, f_H_11] = tfestimate(squeeze(tauIn), squeeze(th1abs), hann(floor(size(squeeze(tauIn), 1) / 10)), [], nfs*4, Fs);
[H_21, f_H_21] = tfestimate(squeeze(tauIn), squeeze(th2abs), hann(floor(size(squeeze(tauIn), 1) / 10)), [], nfs*4, Fs);

% Estimated tf smoothing, by movmean with [x,x] values
H_11 = movmean(H_11,[2,2]);
H_21 = movmean(H_21,[2,2]);

figure
H11sys  = frd(H_11, 2*pi*f_H_11);
bode(H11tf, H11sys, 2*pi*f_H_11), grid on, grid minor
legend('Original H11 tf', 'Identified model', 'location', 'southwest')

figure
H21sys  = frd(H_21, 2*pi*f_H_21);
bode(H21tf, H21sys, 2*pi*f_H_21), grid on, grid minor
legend('Original H21 tf', 'Identified model', 'location', 'southwest')

[Paul]

Hi Sam,

f_H_11 and f_H_21 are in Hz as returned from the calls to tfestimate. What happens if you try

H11sys  = frd(H_11, f_H_11*2*pi);
bode(H11tf, H11sys, f_H_11*2*pi), grid on, grid minor

[Sam Chak]

Hi @Paul

Thank you for your sharp observations. I admit that I seldom use the tfestimate() function from the Signal Processing Toolbox, so I wasn't aware that the unit of frequency is in Hz.

[Paul]

Just to be clear, whether or not the frequency is in Hz depends on how tfestimate is called.

[Paul]

@Sam Chak

Do you have a link to reference for the pole placement technique used here?

[Sam Chak]

My pleasure, @Paul. I am uncertain whether the technique is 100% identical, but they do share certain similarities.

Zítek, P., FiŜer, J. and Vyhlídal, T. (2010) ‘Ultimate-frequency based dominant pole placement’, IFAC Proceedings Volumes, 43(2), pp. 87–92. doi:10.3182/20100607-3-cz-4010.00017.

Answer 1

[Paul]

Hi GIOVANNI,

I don't understand what the math is trying to do or how f and g will be used, so this answer may be limited in scope.

One difference between the two methods is that this line, in Method 1

r = evalfr(H, p(i))*B;

revalutes H at the pole location p(i) in the complex plane. Each p(i) does have a real part.

In Method 2, these lines

H1_interp = interp1(f_H_11, H_11, f_target, 'linear', 'extrap');
H2_interp = interp1(f_H_21, H_21, f_target, 'linear', 'extrap');
r = [H1_interp; H2_interp];  

are evaluating the H_11 and H_21 at a point on imaginary axis, so would be different than r from Method 1.

It looks like the p(i) are very close to the resonances, so maybe small differences in where H is evalated can lead to large differences in r.

[GIOVANNI]

Hi, thank you for the answer

You don't need to understand the math behind those gains, that's only a piece of the project I'm doing.

My doubt is how to make the gains look similar to method 1, so if possible show me how to fix the code to evaluate H in the complex plane like in method 1 but using H_11 and H_21

[Paul]

H_11 and H_21 are defined for values of frequency along the positive, imaginary axis. If you want to evaluate them in the complex plane at locations with non-zero real part, I suppose you can use the data in H_11 and H_21 to try to estimate an actual LTI system, perhaps using tfest or ssest.

Note also that Method 1 evaluates H at four distinct p(i), but Method 2 at only two p(i), because it uses abs.

[GIOVANNI]

The problem is that i'm trying to avoid tfest because it doesn't detect antiresonance close to resonance at 6 Hz, is there any other way to calculate f and g close to 1st method?

[Paul]

tfest seems to work o.k. for H_11 with the following considerations:

Apparently tfest doesn't like the first point in the frequency vector at zero (generates a warning), so get rid of that.

Use a weighting filter to tell tfest that the resonance at ~35 rad/sec is important. See tfestOptions

The high frequency portion of H_11 that increases to a constant value throws off the fit, so get rid of that.

Maybe something similar will work for H_21.

Don't know if this will help for your problem, but mabye it's a start.

% Upload I/O data
%load 'test_identification_analytical_model_.mat'
load 'test_identific...al_model_.mat'
H(1,1)

ans =
 
           113 s^2 + 201.6 s + 1.417e05
  ----------------------------------------------
  s^4 + 3.626 s^3 + 1478 s^2 + 2636 s + 7.984e04
 
Continuous-time transfer function.

% Data measured from simulated model
th1abs = theta1_rad_MEAS.signals.values;   
th2abs = theta2_rad_MEAS.signals.values;
tauIn = tauInput.signals.values;
nfs = 4096*4;   % Number of samples for FFT
Fs = 1000;   % Sampling frequency
% Transfer functions estimate from input 1 (tauIn) to outputs 1 and 2 (th1abs and th2abs)
[H_11, f_H_11] = tfestimate(squeeze(tauIn), squeeze(th1abs), hann(floor(size(squeeze(tauIn), 1) / 10)), [], nfs*4, Fs);
[H_21, f_H_21] = tfestimate(squeeze(tauIn), squeeze(th2abs), hann(floor(size(squeeze(tauIn), 1) / 10)), [], nfs*4, Fs);
% Estimated tf smoothing, by movmean with [x,x] values
H_11 = movmean(H_11,[2,2]);
H_21 = movmean(H_21,[2,2]);
index = f_H_11 > 0 & 2*pi*f_H_11 < 90;  % frequencies that are important
H_11 = frd(H_11(index),2*pi*f_H_11(index));
H_21 = frd(H_21,2*pi*f_H_11);

opt = tfestOptions;
opt.WeightingFilter = 'inv';
sysh11 = tfest(H_11,4,2,opt);
bode(H(1,1),H_11,sysh11,2*pi*f_H_11)

●