Estimating transfer function of a system with input and output values as well as kp and ki values

Question

Keaton Looper 2023 年 11 月 16 日

0
リンク

この質問への直接リンク

https://jp.mathworks.com/matlabcentral/answers/2048252-estimating-transfer-function-of-a-system-with-input-and-output-values-as-well-as-kp-and-ki-values

コメント済み: William Rose 2023 年 12 月 1 日

I am trying to estimate the transfer function of a heat transfer system I am studying. I have the input voltage value and the output corresponding voltage for the temperature sensor. I am using set Kp and Ki values for the system. Is there a way for matlab to calculate the transfer function based and the input, output, and kp and ki values?

1 件のコメント
-1 件の古いコメントを表示-1 件の古いコメントを非表示

Sam Chak 2023 年 11 月 22 日

Hi @Keaton Looper

Based on your description of the Thermal System, I think that the block configuration looks as shown in the image below. To make it easier for @William Rose to guide you in identifying the Thermal System, can you provide the data of Control Voltage in Scope 1 (Vc) and the Output Voltage in Scope 2 (Vo)? I presume that the values for Kp and Ki are fixed.

Terminology:

Tr = Reference temperature
Vr = Reference voltage
Ve = Error voltage
Vc = Control voltage (measured input)
Vo = Output voltage (measured output)
u = Heat from the heating element (energy supplied to the thermal system)
y = Output temperature of the thermal system

サインインしてコメントする。

サインインしてこの質問に回答する。

Answer 1

William Rose 2023 年 11 月 16 日

0
リンク

この回答への直接リンク

https://jp.mathworks.com/matlabcentral/answers/2048252-estimating-transfer-function-of-a-system-with-input-and-output-values-as-well-as-kp-and-ki-values#answer_1354317

@Keaton Looper,

Yes you can.

You can estimate the overall transfer function from Vin(t) to Vout(t) by using tfestimate(), if you have the signal processing toolbox. If you don't have that toolbox, then compute the Fourier transforms of Vin and Vout. Then compute complex ratio Vout(f)/Vin(f) at each frequency. This the (estimate of) the overall transfer function.

Since you mention kp and ki, you must referring to a negative feedback control system. Maybe you want to estimate the open loop transfer function. Let's call it G(f). The overall transfer function which you computed above is

.

You know H(f), right? It is

, therefore

. So you use the known H(f) and the estimated

, and you solve for G(f). That leaves a little for you to work out on your own. Good luck!

35 件のコメント
33 件の古いコメントを表示33 件の古いコメントを非表示

William Rose 2023 年 11 月 17 日

編集済み: William Rose 2023 年 11 月 17 日

[edit: replace figure with a version that displays correctly]

@Paul, thank you, you are right. I did not remember PID controller architecture correctly.

@Keaton Looper, here is a standard diagram of a PID controller, which you can find repeated in many locations.

The setpoint is x(t), the output is y(t). H=1. The process is G, and that is what you want to find, I think. You will first estimate the overall transfer function, T(s)=Y(s)/X(s), as described above. It follows from the diagram above that Error=e(t)=x(t)-y(t). Therefore, in the Laplace domain, E(s)=X(s)-Y(s). You know kp and ki. I assume kd=0. It also follows from the diagram above that

Y(s) = (X(s)-Y(s))*(kp+ki/s)*G(s)

Rearrange to get a formula for the closed loop transfer function, T(s)=Y(s)/X(s).

Then do the algebra to find the formula for G(s) in terms of T(s).

You will first estimate T(s), the closed loop transfer function, from the measured x(t) and y(t). Then you use the formula for G(s) in terms of T(s) (along with your knowledge of kp and ki) to estimate G(s). Of course

.

William Rose 2023 年 11 月 20 日

@Keaton Looper,

from our earlier discussion

T=Y/X closed loop transfer function

From the system architecture: Y=(X-Y)*H*G where H=kp+ki/s

therefore Y(1+GH)=GH*X, therefore T=Y/X=GH/(1+GH)

You determine T by measuring Y and X, or maybe you have T from some other source. YOu want to estimate G. You know H=kp+ki/s=kp+ki*2*pi/(if) where f=frequency in Hz, and i=sqrt(-1).

Therefore T+GH*T=GH, therefore T=GH-GH*T, therefore T=(1-T)*GH, therefore T/(1-T)=GH, therefore

G=T/(H*(1-T))

That's what I get. Check my work, in case I made a mistake.

Let's assume the formula above is correct, and let us use the formula to estimate G.

T should be a vector of complex numbers. Let us assume that it has length N and that the frequencies associated with T are

f=(0:N-1)*df;

and let's assume you have assigned values to kp and ki. Then it is tempting to define

H=kp+2*pi*ki./(f*1i);

The trouble with the above is that it will blow up at f=0. Therefore I recommend skipping 0 frequency:

f1=f(2:end);

T1=T(2:end);

and define

H1=kp+2*pi*ki./(f1*1i);

H1 is now a complex vector of length N-1. Then define

G1=T1./(H1.*(1-T1));

G1 should now be a complex vector of length N-1, and the corresponding frequencies are f1.

I have not checked my suggestions above by actually writing code. Caveat emptor :).

William Rose 2023 年 11 月 22 日

MATLAB Online で開く

tempControllerData.txt

@Keaton Looper,

Here is a plot of your data. I am referring to column 1 as Vc (output from controller) and column 2 as Vo (output from system), in agreement with the diagram provided by @Sam Chak.

data=load('tempControllerData.txt');

fs=1; % sampling rate (Hz)

t=(1:length(data))/fs; % time vector

Vc=data(:,1); Vo=data(:,2);

plot(t,Vc,'r.',t,Vo,'b.');

grid on; xlabel('Time (s)'); ylabel('Temperature');

legend('Vc','Vo');

What is the cause of the disconinuity in Vc at t=12? Did the setpoint change, from something else to 50?

You said you know kp and ki (and kd=0). What are the values?

At the beginning, Vc is positive (Vc~=10), and yet Vo is trending down. Why? I am asking because initial downward trend of Vo suggests (a) that the system is not in a steady state at t=0, or (b) there is a significant unmodelled external disturbance, such as cooling, or (c) the system has inverse response, like the center of mass of a rear wheel drive car, or the water level in some boilers. Another possibility that could explain the trends above is that there is a sign reversal in the system (which is not the same as inverse response). My guess is that the answer is (a) and (b), but your opinion is what matters here.

The data in the figure above make me believe that frequency domain analysis is not going to be useful for this system. I say that because signals above to not have much power at a good range of frequencies, and because te data above represent only a fraction of one transient response. Frequency domain analysis is more likely to give useful results when the measured signals include responses to multiple disturbances or setpoint changes. But maybe, depending on your answers to the previous quesiotns, frequency response will be useful.

ALso, as I alluded to above, the most significant transient in the data is the step down in Vc at t=12. This is assoicated with a rise in Vo. Is the gradual rise in vo, after t=12, due to the decrease in Vc relative to its intitial value, or is the rise in Vo due to the fact that Vc is greater than zero?

More data would be great, if you can get it. For example, record Vc and Vo for a longer period of time, while you change the set point from 40 to 50 to 45 etc, several times, with varying time periods between the changes. And record or take note of the set point changes, of course. But I uderstand you may not be able to get such data.

William Rose 2023 年 11 月 23 日

@Keaton Looper,

Here is a modified version of the diagram posted by @Sam Chak. It shows the two signals you have posted, where vc(t)=Heater and vo(t)=T2. (If the diagram doesn't correctly show the location of the signals, then please explain.)

I think you want to estimate

. Since you have vc(t), which comes after H(s) in the diagram, you do not have to do complex algebra to remove the effect of H(s). Therefore you don't need to know kp or ki, and the estimate of G(s) is simply

where

is the estimate of the cross power spectrum

and

is the estimate of the power spectrum

.(1)

What happens to vc(t) if you turn off the heater (vo=0) for a long time? I am asking because the implicit assumption of the frequency domain method is that the output goes to zero if the input is zero. For many systems, this is not a good assumption, and therefore the DC term in

is ignored.

Compute

with cpsd(vo,vc,...). (If you use cpsd(vc,vo,...), the phase will be off by a factor of -1.) Compute

with cpsd(vc,vc,...) or, equivalently, with pwelch(vc,...).

You might wonder: "Why not just do Ghat=Vo(f)./Vc(f), where Vo(f) and Vc(f) are the one-sided Fourier transforms of vo(t) and vc(t)?" The method "Ghat=Sxy/Sxx" is optimal in a least squares sense, as explained in various textbooks. If you only use a single segment to estimate Sxy and Sxx, then Ghat=Sxy/Sxx and Ghat=Y(f)/X(f) give the same answer. But if you estimate Sxy and Sxx by chopping x(t) and y(t) into separate smaller segments, then the two approaches give different answers. The default behavior of cpsd() and pwelch() is to chop the time domain signal into smaller segments. For your data, I would use a single segment, since your data record is so short. But I still recommend using the approach Ghat=Sxy/Sxx, so that, if you have a longer record in the future, you can easily adapt your code to use the best method.

William Rose 2023 年 11 月 23 日

MATLAB Online で開く

tempControllerData.txt

@Keaton Looper,

data=load('tempControllerData.txt');

vc=data(:,1); % input

vo=data(:,2); % output

fs=1; % sampling rate (Hz)

N=length(vc); % number of samples

t=(1:N)/fs; % time (s)

% compute spectra

[Sco,f]=cpsd(vo,vc,N,N/2,N,fs); % estimate Vc*Vo (* indicates compl.conj)

[Scc,f]=cpsd(vc,vc,N,N/2,N,fs); % estimate Vc*Vc

Ghat=Sco./Scc;

% plot results

figure;

subplot(211), plot(t,vc,'r.',t,vo,'b.')

xlabel('Time (s)'); ylabel('Temp.'); grid on; legend('vc','vo')

subplot(413), plot(f,abs(Ghat),'-k.')

ylabel('|Ghat|'); grid on

subplot(414), plot(f,unwrap(angle(Ghat))*180/pi,'-k.')

xlabel('Frequency (Hz)'); ylabel('phase(Ghat)'); grid on

I do not think the magnitude and phase spectra in the lower two panels of the figure above give an accurate representation of the system G(f), because I think the dominant features of the time domain signals are due to significant unmodeled disturbances and events. For example, |Ghat(f)| has prominent peaks at 0.08, 0.17, 0.25 Hz, etc. These peaks are due to the 12 second long pulse (1/(12 s) = 0.083 Hz) at the beginning of vc(t), which is probably due to an unmodelled change in the setpoint.

I think an approach that is more likely to give useful insight into the system is to model this system in the time domain. That approach will only give useful results if you answer the questions I asked earlier, so that the time domain model can reflect your answers. Those questions include:

Why is vo(t) positive, and trending down, at the start?

What causes the discontinuous jump in vc(t) at t=12?

What would happen to vo(t) if vc(t)=0 for a long time?

William Rose 2023 年 11 月 24 日

MATLAB Online で開く

tempControllerData2.txt

@Keaton Looper,

Plot your data2 (plot appears below):

data=load('tempControllerData2.txt');
N=length(data);
fs1=1;  % experiment sampling rate (Hz) 
t1=(0:N-1)/fs1;  % experiment time vector (s) 
x=data(:,1); T=data(:,2);
plot(t1,T,'bo',t1,x,'ro'); grid on; hold on
xlabel('Time (s)'); ylabel('Temperature');
legend('Texp','Heater')

Do you have a guess for the transfer function? I'm not great at transfer functions, so I would rather write a simple differential equation, by thinking about how heaters work and how objects or tanks of liquid gain and lose heat. The differential equation below is a first guess:

which says the net rate of change of temperature is the sum of heat heat loss to the environment plus heat gain from the heater. x is the input to the heater,

is a time constant for heat exchange with the environment, and

is a time constant for the heater. Note that the system above trends toward T=0 as T goes to infinity, if the heat is off (x=0). That makes sense if the ambient temperature is zero.

The system above has the Laplace-transform representation

therefore the transfer function from X to T is

You can estimate

and

from your data with lsqnonlin(), or by some other fitting process. Note that the system above approaches T(t)/X(t)=tau1/tau2, as t goes to infinity. The experimental data you provided (data2) shows that the sytem is prettty close to steady state at t=85 s, with T=48 and x=3.4. Therefore I choose to define tau2=tau1*(3.4/48), for initial testing of this model. Then there is just one parameter to estimate: tau1. When you play with the model above, you will find that it does not do the job. For example, with tau1=90 s, you get tau2=6.38 s, and this:

fs2=10; % simulation sampling rate (Hz)

tau1=90; % time constant 1 (try different values)

tau2=tau1*3.4/48; % time constant 2 (see text for explanation)

t2=(0:N*fs2/fs1-1)/fs2; % simulation time vector (s)

x2=interp1(t1,x,t2); % heater signal, interpolated to the simulation sampling rate

Tsim=Tsim1(t2,x2,tau1,tau2,T(1));

plot(t2,Tsim,'b-',t2,x2,'r-');

legend('Texp','Heater','Tsim','HeatSim','Location','southeast')

function T=Tsim1(t,h,tau1,tau2,T0)

% Compute simulated temperature (Model 1)

% Inputs

% t=time vector (s), h=heater signal

% tau1, tau2=time constants for heat echange with environment and for heater

% T0=initial temperature

% Output

% T=vector of simulated temperatures

N=length(t);

T=T0*ones(N,1);

dt=(t(N)-t(1))/(N-1);

for i=2:N

T(i)=T(i-1)+dt*(-T(i-1)/tau1+h(i-1)/tau2);

end

The plot above shows that we got the beginning and the end right, but not the dynamics. Try different values for tau1 in the code above, to get a feel for the system.

The result above convinces me that we need another term in the transfer function,or (equivalently) another piece in the differential equation model.

William Rose 2023 年 11 月 24 日

@Keaton Looper,

I said earlier that I think we need another term in the transfer function, or another piece in the differential equation model, to model your data2. The model that did not work well was

which corresponds to

and which has transfer function

I recommend we low-pass x(t) to get a better fit. Physically, we can think of this as a 2-stage heater. In stage 1, x(t) heats an intermediate heating fluid, or outer jacket, which has temperature Ti(t). In stage 2, the intermediate fluid heat the final tank of fluid, whose temperature is T(t). Now the differential equations are

Stage 1:

and Stage 2:

which can be combined to make

whose transfer function is

So now we have a pair of first order differential equation (stage 1 and stage 2 equations, above) and we have the corresponding transfer function. The next step is to find values of the unknown constants tau_x, tau1, and tau2. The transfer function has DC gain=tau1/tau_x, which means in the steady state, T/x=tau_1/taux. Therefore we define taux=tau_1*xfinal/Tfinal=tau_1*(3.4/48). This leaves us with two unknowns, tau1 and tau2. Write a new simulation function, Tsim2(...), that is a modified version of Tsim1(...), and use lsqcurvefit() to fit tau1 and tau2.

William Rose 2023 年 11 月 27 日

編集済み: William Rose 2023 年 11 月 28 日

MATLAB Online で開く

@Keaton Looper,

[edit fix typos] [edit: Fix description of the inputs to Tsim2(), in my text. The code is not changed.]

I suggested that Tsim2 would be like Tsim1, but with a lowpass filter. Wit that change, the differential equations for Tsim2, and the transfer function for Tsim2, were given in a comment I made above.

A Matlab function that implements Tsim2 is below.

function T=Tsim2(tau,texp)
% Compute simulated temperature (Model 2)
% Inputs
%   tau=[tau1,tau2,taux]=time constants (s)
%   texp=time vector for measured data (s)
% Output
%   T=vector of simulated temperatures at times texp
% Global variables used
%   xexp=heater signal
%   Texp0=initial measured temperature
global xexp Texp0
tau1=tau(1); tau2=tau(2); taux=tau(3);
dt=0.05;                            % dt for simulation
tsim=texp(1):dt:texp(end);          % simulation time vector
Nsim=length(tsim);                  % number of simulated points
Ti  =Texp0*ones(Nsim,1);            % allocate array for Ti
Tsim=Texp0*ones(Nsim,1);            % allocate array for T
xsim=interp1(texp,xexp,tsim);       % x at times tsim
for i=2:Nsim
    dTidt  = xsim(i-1)/taux-Ti(i-1)/tau1;
    dTdt   = (Ti(i-1)-Tsim(i-1))/tau2;
    Ti(i)  = Ti(i-1)  +dt*dTidt;    % Ti, simulated
    Tsim(i)= Tsim(i-1)+dt*dTdt;     % T, simulated
end
T=interp1(tsim,Tsim,texp)';         % Tsim at times texp
end

Function Tsim2(tau, texp) returns a vector which is the s.imulated system output - the temperature. The input tau is a vector with three components: tau1,tau2,taux. The input texp is a vector of times at which we want to know the simulated temperature. These times are the times when temp[erature was measured. This will allow us to directly compare the measured temps to the simulated temps. The output from Tsim2() is a vector with the same size as the experimentazlly measured temps. I set it up this way because this is what the function lsqcurvefit() wants. lsqcurvefit() will call Tsim2() many times, adjusting the three elements of tau each time, until it finds a set of 3 taus that minimize the error between the simulated and mesaured temperatures

Tsim2() uses the global variables x=input waveform and Texp0=initial measured temperature. I realize that global variables in functions are not recommended, because of the performance penalty, but for this simple system, ther global variables in the function are not a problem.

My Tsim2() is somewhat more complicated than might be expected. The reason it is complicated is that your experimental data is sampled at dt=1 s intervals, and I am not comfortable doing integration by the Euler method for this system with such large time steps. In other words, I suspect the Euler method, for this system, will have significant inaccuracy if implemented with dt=1 s. Possible alternatives include: 1. Don't use Euler, use ode45() instead. I avoided this, because I did not want to have to explain ode45(). 2. Interpolate the input to a shorter sampling interval (I chose dt=0.05 s), run the Euler method at the shorter sampling interval, then interpolate the results to the time points specified in texp. I chose option 2.

Keaton Looper 2023 年 11 月 30 日

i am tyring to graph the data and the estimates but i am unable to get it to graph based on the new Tsim2 function

William Rose 2023 年 12 月 1 日

Keaton, can you share what you tried that didn’t work? Thanks. I am traveling and don’t have a computer with Matlab.

サインインしてコメントする。

Answer 2

William Rose 2023 年 11 月 24 日

0
リンク

この回答への直接リンク

https://jp.mathworks.com/matlabcentral/answers/2048252-estimating-transfer-function-of-a-system-with-input-and-output-values-as-well-as-kp-and-ki-values#answer_1358827

Keaton, what do you get if you analyze the new data with code used to analyze the earlier data? What differences between the transfer functions of the old and new data stand out to you?