coefficient determination of two variables

Hello everyone,

I would like to ask you some guidelines about how to solve this equation in matlab. This is an equation for an admittance of the open ended coaxial probe.

I would like to find the coefficients $\alpha$ and $\beta$ in the equation based on my measured data, in which those parameters are known: $\epsilon_r$, s, a and Y. The upper bound (N,M,P,Q) we can choose. For example, N=M = 4; P=Q=8.

I think some kind of fitting function might address my question, but I am really new to Matlab and have no ideas to do that.

Thank you very much,

Greetings

Answer 1

[Matt J]

lsqcurvefit would be one choice,

https://www.mathworks.com/help/optim/ug/lsqcurvefit.html

Note that your equation can be reorganized as a linear equation by multiplying through by the denominator expression involving the betas. That would likely be a good way fo developing an initial guess for the parameters.

[Duy Hai Nguyen]

Dear Matt,

thank you for your suggestions. I try to write first codes based on this, but it is still not successful yet. Do you have any ideas what should I change? Data file I also included.

% define upper bound of the sum series
N = 3; P = 2; M = 2; Q = 2;

% import measured data
importLoad = importdata('26062020_IndiThera_gufVNA_matlab_calKit3m5_water_outpdata_1.txt');
freq = importLoad.data(:,1)*(1e9); %% frequency term, equivalent to 's' in the eqn, unit: Hz
Ydata = complex(importLoad.data(:,2),importLoad.data(:,3)); %% complex admittance Y

% define data set for epsilon_r (eps_ref)
a = 0.25*(1e-3);
eps_1 = 78.32;
eps_inf = 4.57;
tau = 8.38e-12;
eps_ref= (eps_1 - eps_inf)./(1+1i*2*pi*freq*tau) + eps_inf;

% define numerator of the Y_model, without the coefficients
%% coefficients will be multiplied in the final function
num = zeros(length(freq),1);
for i=1:length(freq)
    for n=1:N
        for p=1:P
            num(i) = num(i) + sqrt(eps_ref(i))^p+(freq(i)*a)^n;
        end
    end
end

% define denominator of the Y_model, without the coefficients
den = zeros(length(freq),1);
for i=1:length(freq)
    for m=1:M
        for q=1:Q
            den(i) = den(i) + sqrt(eps_ref(i))^p+(freq(i)*a)^n;
        end
    end
end

% define the Y_model equation
%% parameters to determine: alpha, beta
%% parameters which are known: freq, eps_Ref
Ymodel = @(alpha, beta,freq, eps_ref) alpha(:,1).*num(:,1)./(1+beta(:,1).*den(:,1)); 

% define initial values of the coefficients
alpha_0 = zeros(length(freq));
beta_0 = zeros(length(freq));

% find the coefficients
cal = lsqcurvefit(Ymodel,alpha_0,beta_0,freq,eps_ref,Ydata);

Thank you very much,

Greetings,

[Matt J]

Ymodel = @(alphabeta,~) alphabeta(:,1).*num./(1+alphabeta(:,2).*den);

% define initial values of the coefficients
alpha_0 = zeros(length(freq),1);
beta_0 = zeros(length(freq),1);

% find the coefficients
alphabeta = lsqcurvefit(Ymodel,[alpha_0,beta_0],[],Ydata);

[Duy Hai Nguyen]

Dear Matt,

thank you very much for your time and fast reply. It worked flawlessly. Unfortunately I am not quite clear about your code. Could you explain the meaning of the following lines:

Ymodel = @(alphabeta,~) alphabeta(:,1).*num./(1+alphabeta(:,2).*den);

What does '~' mean in this case? Why don't we use 2 separated array of alpha and beta instead?

alphabeta = lsqcurvefit(Ymodel,[alpha_0,beta_0],[],Ydata);

What does '[ ]' mean? How can the fitting function understand which variables he will use?

When I set the initial guess of alpha and beta like above, beta result in all zeros. When I change initial values to 1 for both coefficients, it results in different values? I thought that we should have unique coefficients, isn't it?

Thank you.

Greetings

[Matt J]

What does '~' mean in this case?

The tilde can be used to tell Matlab that a function input argument or output argument will not be used,

https://www.mathworks.com/help/matlab/matlab_prog/ignore-function-inputs.html

Since you are not using a 2nd input anywhere in your definition of Ymodel, it makes sense here.

Why don't we use 2 separated array of alpha and beta instead?

Because the documentation for lsqcurvefit tells you that they should not be separated.

"fun is a function that takes two inputs: a vector or matrix x (the unknowns), and a vector or matrix xdata."

In other words, lsqcurvefit expects all unknowns to be part of one single array, both as input to Ymodel and in the initial guess. It also expects all of your xdata to be part of a single array as well.

What does '[ ]' mean?

It is an empty matrix. You could have put anything there. Because of the '~', it will be passed to Ymodel as the 2nd input argument and ignored.

I thought that we should have unique coefficients, isn't it?

I see no reason why the solution would be unique. There is no convexity in your least squares cost function that I can see. However, I did give you a recommendation in my initial post above how you could generate a more educated initial guess - by reorganizing the Ymodel equation in linear form and solving.

[Duy Hai Nguyen]

Dear Matt,

thank you for the explanation. I will try to follow your suggestion. Another question, are there any "smart" way to write the sum in the function handle instead of making two for loops for the numerators and denominator as I did. I tried several ways but still have no ideas how to do that.

Thank you very much.

Greetings,

[Matt J]

Yes, you can eliminate many of the internal loops by rewritting the numerator and denominator expressions in Ymodel as matrix quadratic forms x.'*A*y.

Ymodel = @(alphabeta,~) modelfunc(alphabeta,freq,eps_ref,a,N,P,M,Q);

function out = modelfunc(alphabeta,freq,eps_ref,a,N,P,M,Q)

     Alpha=reshape(alphabeta(:,1),N,P));    %unpack into matrix form
     Beta=reshape(alphabeta(:,2),M,Q+1));
    
     Ns=numel(freq);
     out=nan(1,Ns);
     
     for i=1:Ns

         sa=freq(i)*a;
         sqrt_epsr=sqrt(eps_ref(i));
         
         san=(sa).^(1:N);
         sam=(sa).^(1:M);
        
         erp=sqrt_epsr.^(1:P);
         erq=sqrt_epsr.^(0:Q);
        
         out(i)=san*Alpha*erp(:)./(1+ sam*Beta*erq(:));
         
     end

 
end

[Duy Hai Nguyen]

Dear Matt,

many thanks.

[Duy Hai Nguyen]

Dear Matt,

sorry to bother you again. When I try your code alone, there is no error. However, when I put the lsqcurvefit to find the alpha, beta. I got error of the size of the matrix. Please see below:

Here is the code for the alphabeta, I guess there should not any changes:

AlphaBeta = lsqcurvefit(y_model,[alpha_0,beta_0],[],y_ref);

And here is the error:

Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that
dimension.

Error in openEndedCoax>modelfunc (line 61)
     Alpha=reshape(alphabeta(:,1),N,P);    %unpack into matrix form

Error in openEndedCoax>@(alphabeta,~)modelfunc(alphabeta,freq,eps_ref,a,N,P,M,Q) (line 49)
y_model = @(alphabeta,~) modelfunc(alphabeta,freq,eps_ref,a,N,P,M,Q);

Error in lsqcurvefit (line 222)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in openEndedCoax (line 58)
AlphaBeta = lsqcurvefit(y_model,[alpha_0,beta_0],[],y_ref);

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

Do you have any ideas?

Greetings,

[Matt J]

Yes, the N and P you've passed in do not agree with the size of alphabeta. You should get in the habit of using dbstop

https://www.mathworks.com/help/matlab/ref/dbstop.html

to probe such errors.