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

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?

[Sam Chak]

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]

@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!

[William Rose]

@Keaton Looper, remember that G(f) and H(f) are complex vectors, in the formulas above.

[Paul]

I think it's much more common for the PI control to be in the forward path, in which case the closed loop transfer function would be H(f)*G(f)/(1 + H(f)*G(f))

[William Rose]

[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 .

[Keaton Looper]

I am using the way the is described above and I am getting a 129x129 complex double of all 0 with now transfer function shown.

[William Rose]

@Keaton Looper, that does not sound good.

Are you doing "dot-divide", i.e. element-by-element divide, of two vectors, at the appropriate points? More later, no time right now.

[Keaton Looper]

What would that look like in matlab code for example? I found the function for G(s) in terms of T(s).

[William Rose]

@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 :).

[Keaton Looper]

So what would I do if I have a set point of 50 constant for the system and a sample rate of 1 samples per second?

[William Rose]

@Keaton Looper,

Do you have the input and output data? If so, post it.

Your original post was "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?"

I assumed based on your quesiton that you would have a recording of the input and output volatges and that they would show some variability with time. If this is not true (i.e. if the input never changes) then you will not have much luck, because there is no power in the input at any frequency, except at f=0.

[Keaton Looper]

The input does change because of the negative feedback loop but there is a set point that the system is trying to reach. I will post the data in a bit.

[Keaton Looper]

is there an overall transfer function to be calculated that shows the numerator and denominator? when i caluclate it it gives me a 128x1 complex double. What do these numbers mean?

[William Rose]

Keaton, I assume you have experimental data and a known kp and ki. I assume you want to estimate the transfer function of the plant from your data. Please post the data you have, and a clear description and diagram showing where in the system the data is from. I can’t answer your specific question about the complex vector yet, because I don’t know what data you used or what formula you used.

[Keaton Looper]

As commented by @Sam Chak the heater has a set poont voltage which is being inserted into the system. The kp and ki values are set constant and the temperature is being read as the heater voltage changes to try to get T2 to the temperture set point.

The data below in orange represent the data measured from the heater after the controller, and the grey represents the data for the temperature T2 being measured as it approaches the set point of 50. The columns are continuous from left to right and have been put side by side so that it could fit in one image.

[Keaton Looper]

Here is the table with the data

[William Rose]

@Keaton Looper, thanks. More later.

[William Rose]

@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.

[Keaton Looper]

I am unable to collect more data but I can post the data from the second trial using the same kp ki and set point conditions.

[William Rose]

@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,...).

1. 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]

@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?

[Keaton Looper]

I can upload data from the trial of the same parameters to see if it has the same discontinuity.

[Sam Chak]

Hi @Keaton Looper

If possible inject a sinusoidal signal at the reference.

t = linspace(0, 1000, 10001);
T = 5*sin(2*pi/1000*t) + 45;    % Sinusoidal time-varying setpoint
plot(t, T), grid on
xlabel('t / seconds')
ylabel('T / °C')

[Keaton Looper]

Here is the second trial of the same testing parameters. This data doesnt have a discontinuity therefore am i able to estimate the overall transfer function for the system?

[Keaton Looper]

Is there still a way for me to get the numerator and denominator of the transfer function representing the system?

[William Rose]

@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
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.

[Keaton Looper]

Maybe it’s missing the part of the system which is the kp and ki value or is that not included?

[William Rose]

@Keaton Looper,

In case it is not obvious, I wrote function Tsim1(...). It returns a vector of simulated temperatures. It uses the differential equations above, integrated by Euler's method. The variables passed to the funciton are described in the comments. I used a time vector with 0.1 s spacing, rather than the 1 s spacing of the experiment, in order to have smaller errors with Euler's method. Since I am only using the output from Tsim1 for visual comparison to the experimental data, it is OK that the experiment and the simulaiton have different time spacing. If you were going to use the output from Tsim1() to compute the squared error between experiment and model, then you'd want the times in Tsim1 to match the times of the experimental data. But then you'd have to consider whether Euler's method for integration is accurate enough at the experimental sampling rate.

[William Rose]

@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.

[Keaton Looper]

I am not familiar with lsqcurvefit() how is it formatted to find tau1 and tau2?

[William Rose]

@Keaton Looper,

Please read the help for lsqcurvefit() to learn how to use it. That's what I do. There are good examples there. Then, if you get stuck, post what you have, describe the error that happens, and I bet someone on this site will be able to help you.

Good luck.

[Sam Chak]

What are the values for Kp and Ki used in the second trial (or Data 2)?

[Keaton Looper]

I think kp was 5 and ki was 0.05

[Keaton Looper]

@William Rose I have having trouble determining what the function would be for Tsim2. Wouldn't it be the same as Tsim1 since its a function of tau1 and tau2?

[William Rose]

@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]

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]

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]

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?

Answer 3

[Branko]

Hi @William Rose

I am also a student and would like to simulate a heating system with a heat pump.

The heat pump is controlled as on/off or by adjusting the inflow curve from -5 C to 5 C.

Can you please help me with instructions on how to start building the system.

Thank you.

Regards,

Branko