pole placement controller gain trouble with identified receptance

Question

GIOVANNI 2025-4-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175993-pole-placement-controller-gain-trouble-with-identified-receptance

评论： Sam Chak 2025-4-5

test_identification_analytical_model_.mat

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

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

f = 1×2

1.0e-06 * -0.3688 -0.3196

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

% 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

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

f = 1×2

0 0

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

Sam Chak 2025-4-5

编辑：Sam Chak 2025-4-5

在 MATLAB Online 中打开

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 2025-4-5

@Sam Chak

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

Sam Chak 2025-4-5

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 2025-4-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175993-pole-placement-controller-gain-trouble-with-identified-receptance#answer_1563171

在 MATLAB Online 中打开

test_identification_analytical_model_.mat

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

GIOVANNI 2025-4-5

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 2025-4-5

在 MATLAB Online 中打开

test_identification_analytical_model_.mat

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)

●

请先登录，再进行评论。

pole placement controller gain trouble with identified receptance

8 个评论
显示 6更早的评论隐藏 6更早的评论

回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

pole placement controller gain trouble with identified receptance

8 个评论 显示 6更早的评论隐藏 6更早的评论

回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论