- 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

**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

63 views (last 30 days)

Show older comments

##### 1 Comment

Sam Chak
on 22 Nov 2023

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:

### Answers (3)

William Rose
on 16 Nov 2023

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 Comments

William Rose
on 17 Nov 2023

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

Paul
on 17 Nov 2023

William Rose
on 17 Nov 2023

Edited: William Rose
on 17 Nov 2023

[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
on 20 Nov 2023

William Rose
on 20 Nov 2023

@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
on 20 Nov 2023

William Rose
on 20 Nov 2023

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
on 20 Nov 2023

William Rose
on 21 Nov 2023

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
on 21 Nov 2023

Keaton Looper
on 21 Nov 2023

William Rose
on 22 Nov 2023

Keaton Looper
on 22 Nov 2023

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.

William Rose
on 22 Nov 2023

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
on 22 Nov 2023

William Rose
on 23 Nov 2023

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
on 23 Nov 2023

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
on 23 Nov 2023

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

Sam Chak
on 23 Nov 2023

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
on 23 Nov 2023

Keaton Looper
on 24 Nov 2023

William Rose
on 24 Nov 2023

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
on 24 Nov 2023

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

William Rose
on 24 Nov 2023

William Rose
on 24 Nov 2023

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
on 24 Nov 2023

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

William Rose
on 25 Nov 2023

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.

Keaton Looper
on 27 Nov 2023

William Rose
on 27 Nov 2023

Edited: William Rose
on 28 Nov 2023

[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
on 30 Nov 2023

William Rose
on 1 Dec 2023

William Rose
on 24 Nov 2023

Branko
on 24 Nov 2023

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

### See Also

### Categories

### Tags

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)