How do I run the function for each value in array?

Hello,

the following code calculates two time constants to correct thermocouple measurements. At the moment it works for one pair of measured temperature at the time x, Tth1 and Tth2. The goal is to import columns of data pairs out of excel and let the code run for each data pair and give the time constants for each pair to get the mean time constant afterwards with least squares method.

Problem is that the code can not handle imported column vectors and just gives an error message:

Error using fmincon
Supplied objective function must return a scalar value.

Error in versuch (line 53)
[tau,fval] = fmincon(fun,tau0,[],[],[],[],[0,0],[],[]);

Everything I tried until now will not give me any results, it just works if I type in each pair seperatly.

An example of the excel data pairs is attached.

clear

global Tth1 Tth2 A B
data= 
%Tth1=1. thermocouple bead temperature
Tth1= 1747.15;
%Tth2=2. thermocouple bead temperature
Tth2= 1704.15;
%d1= 1. bead diameter (m) 
d1=0.000125;
%d2= 2. bead diameter (m)
d2=0.000325;
%epsilon= bead emissivity(assumed 0.95 for soot coated)
epsilon=0.95;
%sigma= Stefan Boltzman constant (W/(m^2*K^4)
sigma=5.67*10^-8;
%rho= bead density(kg/m^2)
rho=21460;
%cp= [J/(kg*K)]
cp=133;

A=(6*epsilon*sigma)/(rho*cp*d1);
B=(6*epsilon*sigma)/(rho*cp*d2);



% calculating the results
fun = @(tau)(Tth1 + tau(1)*A*Tth1^4 - (Tth2 + tau(2)*B*Tth2^4))^2;
tau0 = [0,0];
[tau,fval] = fmincon(fun,tau0,[],[],[],[],[0,0],[],[]);

% display of results
disp(['tau1 = ',num2str(tau(1))])
disp(['tau2 = ',num2str(tau(2))])

 
% check
% y1 and y2 should be identical at best
y1 = Tth1 + tau(1)*A*Tth1^4;
y2 = Tth2 + tau(2)*B*Tth2^4;
deviation = abs(y1-y2);
disp(['Abweichung: ',num2str(deviation)]);

deviation

[Matheo Schaaf]

Thanks for the help! I'm relatively new to MATLAB and for now I'm unsure how to write the drivers script so that each of my 34000 values gets used.

Answer 1

[dpb]

clear
global Tth1 Tth2 A B
data= 
...

Therein lies your problem -- you wrote a script instead of a function -- get rid of the clear statement; wrap the code in a function with the needed inputs passed as arguments.

Then, write a driver script that calls the function for each pair of inputs in turn, saving the results returned as outputs of the function, NOT as global variables.

You undoubtedly could also vectorize the function itself to accept an array or pair of vectors and do the looping internally or using MATLAB builtin vectorized operations, but I didn't look at the code enough to see how much trouble that might turn out to be; the first option is almost a trivial exercise.

[dpb]

"... I'm unsure how to write the drivers script so that each of my 34000 values gets used. "

Read the data from wherever they are into an array and then loop over them...

tData=readmatrix('yourDataFile.dat');   % import the data
N=size(tData,1);                        % the number of observations
tau(N,3)=0;                             % preallocate results (2 tau, one difference)
for i=1:N
  tau(i,:)=yourFunction(tData(i,:));    % call your function with the two T values
end

your function as above will return a 3-vector of the two tau values and difference from each pair of

inputs; modify it and the return as see fit to match your desires...

function res=chooseYourFunctionName(T)
    % assign input T vector into local variables
    
    % dpb changes to original code...
    Tth1=T(1);
    Tth2=T(2);
    % -- to here.  That's it other than removing globals and the clearstuff
    
    %d1= 1. bead diameter (m) 
    d1=0.000125;
    %d2= 2. bead diameter (m)
    d2=0.000325;
    %epsilon= bead emissivity(assumed 0.95 for soot coated)
    epsilon=0.95;
    %sigma= Stefan Boltzman constant (W/(m^2*K^4)
    sigma=5.67*10^-8;
    %rho= bead density(kg/m^2)
    rho=21460;
    %cp= [J/(kg*K)]
    cp=133;
    A=(6*epsilon*sigma)/(rho*cp*d1);
    B=(6*epsilon*sigma)/(rho*cp*d2);
    % calculating the results
    fun = @(tau)(Tth1 + tau(1)*A*Tth1^4 - (Tth2 + tau(2)*B*Tth2^4))^2;
    tau0 = [0,0];
    [tau,fval] = fmincon(fun,tau0,[],[],[],[],[0,0],[],[]);
    
    % display of results
    % won't want several thousand rolling off screen in batch mode...
    % could have an optional debug flag variable to control...
    %disp(['tau1 = ',num2str(tau(1))])
    %disp(['tau2 = ',num2str(tau(2))])
    
    
    % check
    % y1 and y2 should be identical at best
    y1 = Tth1 + tau(1)*A*Tth1^4;
    y2 = Tth2 + tau(2)*B*Tth2^4;
    deviation = abs(y1-y2);
    %disp(['Abweichung: ',num2str(deviation)]);
    
    % package output for return to caller...
    res=[tau deviation];
end

Above only passes in the T array and assumes all constants are, indeed, constant. If there's reason to change some of those routinely as well, then add those to the argument list as (perhaps optional) arguments.

Read up in the section of programming about functions -- it'll explain all of this and more with examples.

[Walter Roberson]

See also arrayfun()

[dpb]

<Is-arrayfun-faster-than-for-loop?> (in which thread you also participated, Walter :) ) shows that depending on the function and whether in m-file or command line, the overhead with arrayfun may be sizable. Probably even with OPs size of some 35K elements the time won't be but fraction of to maybe seconds (depending on how the fmincon handles the solution) so I opted for the loop...besides being less abstract for the newcomer.

I thought about but hadn't tried to see if perhaps fsolve wouldn't perhaps be faster.