how can i estimate parameters that must be real from a complex model?

Hi, I have a problem with estimating the coefficients of a complex nonlinear model. In practice the model is

Measurements were carried out, as a function of the frequency f, obtaining εv, which is composed of the real part and the imaginary part. The goal is to estimate the parameters εr, vfw, vb, sigma that fit the non-linear model. I have used several functions like lsq or lsqnonlin. My problem is to estimate these parameters, which must be real, at the same time for the real and the imaginary part. I've read that I have to set an objective function that takes into account the real part and the imaginary part, but I can't figure it out. Can anyone kindly help me

Answer 1

[Torsten]

f = [...];     %your data row vector;
eps_v = [...]; %your data row vector;
%Initial values for parameters to be fitted
eps_r0 = ...;
vfw0 = ...;
vb0 = ...;
sigma0 = ...;
p0 = [eps_r0,vfw0,vb0,sigma0];
% Call nonlinear optimizer
p = lsqnonlin(@(p)fun(p,f,eps_v),p0)
    
function res = fun(p,f,eps_v)
  eps_r = p(1);
  vfw = p(2);
  vb = p(3);
  sigma  =p(4);
  eps_v_model = eps_r + vfw*(4.9+75./(1+1i*f/18)-1i*18*sigma./f)+...
                vb*(2.9+55./(1+(1i*f/18).^(0.5)));
  res1 = real(eps_v_model-eps_v);
  res2 = imag(eps_v_model-eps_v);
  res = [res1,res2];
end  

[daniele gallo]

thank you for your answer, i just have some doubts to ask you, you have to excuse me but i'm new with matlab. The measurement data are contained in a csv file, which is composed as follows: first frequency column, second real eps_v column, third imaginary eps_v column. is it correct to load this way? f = frequencies

epr = real part;

epm = - (partImag);

eps_v = [epr epm];

also if I wanted to use the lsqcurfit function, would it change something in the setting you made above?

[Torsten]

If you already have separated eps_v in real and imaginary part, then hold them as they are.

    f = [...]; %your data row vector;
    epr = [...]; %real part of your data row vector;
    epi = [...]; %imaginary part of your data row vector;
    eps_v = [epr,epi];
    %Initial values for parameters to be fitted
    eps_r0 = ...;
    vfw0 = ...;
    vb0 = ...;
    sigma0 = ...;
    p0 = [eps_r0,vfw0,vb0,sigma0];
    % Call nonlinear optimizer
    p = lsqcurvefit(@fun,p0,f,eps_v)
    
    function res = fun(p,f)
      eps_r = p(1);
      vfw = p(2);
      vb = p(3);
      sigma  =p(4);
      eps_v_model = eps_r + vfw*(4.9+75./(1+1i*f/18)-1i*18*sigma./f)+...
                    vb*(2.9+55./(1+(1i*f/18).^(0.5)));
      res = [real(eps_v_model),imag(eps_v_model)];
    end  

[daniele gallo]

I tried to run the code, but the function calculates parameters that do not fit the model well, they are all null. I load the code and the plots obtained,I have also attached the data if you want to try to check by yourself. can you help me understand where am i wrong? you are very nice

f = frequenzee; %your data row vector;
    epr =partReale %real part of your data row vector;
    epi = -(parteImag); %imaginary part of your data row vector;
    eps_v = [epr,epi];
    %Initial values for parameters to be fitted
  % epv=epr-epi.*1i;
   eps_r0 = 0;
    vfw0 = 0;
    vb0 = 0;
    sigma0 = 0;
   options = optimset('Display','iter','FunValCheck','on', ...
                   'MaxFunEvals',Inf,'MaxIter',Inf, ...
                   'TolFun',1e-6,'TolX',1e-6);
   paramslb = [0 0 0 0];  % lower bound
paramsub = [ ];

    p0 = [eps_r0,vfw0,vb0,sigma0];
    % Call nonlinear optimizer
    p = lsqcurvefit(@fun,p0,f,eps_v,paramslb,paramsub, options)
    stima= p(1) + p(2).*(4.9+75./(1+1i*f/18)-1i*18.*p(4)./f)+...
                   p(3).*(2.9+55./(1+(1i*f/18).^(0.5)));
    plot(f,real(stima))
    hold on
    plot(f,epr)
    figure
    plot(f,imag(stima))
    hold on
    plot(f,epi)
    function res = fun(p,f)
      eps_r = p(1);
      vfw = p(2);
      vb = p(3);
      sigma  =p(4);
      eps_v_model = eps_r + vfw.*(4.9+75./(1+1i*f/18)-1i*18.*sigma./f)+...
                    vb.*(2.9+55./(1+(1i*f/18).^(0.5)));
      res = [real(eps_v_model),imag(eps_v_model)];
    end  

[Torsten]

Taking

epi = parteImag.' 

instead of

epi = -parteImag.'

works. Maybe you have to reconsider whether the latter is correct despite your setting.

load frequenzee.mat
load parteImag.mat
load partReale.mat
f = frequenzee.';
epr = partReale.' %real part of your data row vector;
epi = parteImag.'; %imaginary part of your data row vector;
eps_v = [epr,epi]
%Initial values for parameters to be fitted
eps_r0 = 1;
vfw0 = 1;
vb0 = 1;
sigma0 = 1;
%options = optimset('Display','iter','FunValCheck','on', ...
%               'MaxFunEvals',Inf,'MaxIter',Inf, ...
%               'TolFun',1e-6,'TolX',1e-6);
paramslb = [0 0 0 0];  % lower bound
paramsub = [Inf Inf Inf Inf ];
p0 = [eps_r0,vfw0,vb0,sigma0];
% Call nonlinear optimizer
p = lsqcurvefit(@fun,p0,f,eps_v,paramslb,paramsub);%, options)
stima= p(1) + p(2).*(4.9+75./(1+1i*f/18)-1i*18.*p(4)./f)+...
       p(3).*(2.9+55./(1+(1i*f/18).^(0.5)));
plot(f,real(stima))
hold on
plot(f,epr)
figure
plot(f,imag(stima))
hold on
plot(f,epi)
end
function res = fun(p,f)
  eps_r = p(1);
  vfw = p(2);
  vb = p(3);
  sigma  =p(4);
  eps_v_model = eps_r + vfw.*(4.9+75./(1+1i*f/18)-1i*18.*sigma./f)+...
                vb.*(2.9+55./(1+(1i*f/18).^(0.5)));  
  res = [real(eps_v_model),imag(eps_v_model)];
end  

[daniele gallo]

thank you for your wonderful posting.

how can i make sure that the estimated parameters respect constraints? for example they must all be positive, in particular it must be that vwf + vb = 1 and 1.5 <eps-r <4.5

[Torsten]

This reduces your model to a model with 3 parameters.

Here is my suggestion:

load frequenzee.mat
load parteImag.mat
load partReale.mat
f = frequenzee.';
epr = partReale.' %real part of your data row vector;
epi = parteImag.'; %imaginary part of your data row vector;
eps_v = [epr,epi]
%Initial values for parameters to be fitted
eps_r0 = 1;
vfw0 = 0.5;
sigma0 = 1;
%options = optimset('Display','iter','FunValCheck','on', ...
%               'MaxFunEvals',Inf,'MaxIter',Inf, ...
%               'TolFun',1e-6,'TolX',1e-6);
paramslb = [1.5 0 0];  % lower bound
paramsub = [4.5 2*pi Inf ];
p0 = [eps_r0,asin(sqrt(vfw0)),sigma0];
% Call nonlinear optimizer
p = lsqcurvefit(@fun,p0,f,eps_v,paramslb,paramsub);%, options)
eps_r = p(1)
vfw = (sin(p(2)))^2
vb = (cos(p(2)))^2
sigma = p(3)
stima= eps_r + vfw.*(4.9+75./(1+1i*f/18)-1i*18.*sigma./f)+...
       vb.*(2.9+55./(1+(1i*f/18).^(0.5)));
plot(f,real(stima))
hold on
plot(f,epr)
figure
plot(f,imag(stima))
hold on
plot(f,epi)

function res = fun(p,f)
  eps_r = p(1);
  vfw = (sin(p(2)))^2;
  vb = (cos(p(2)))^2;
  sigma = p(3);
  eps_v_model = eps_r + vfw.*(4.9+75./(1+1i*f/18)-1i*18.*sigma./f)+...
                vb.*(2.9+55./(1+(1i*f/18).^(0.5)));  
  res = [real(eps_v_model),imag(eps_v_model)];
end  

[daniele gallo]

if I understand correctly do you want to exploit the fundamental goniometric relationship?

having no useful information, given the experimental nature, as an initial value of the parameters p0 = [0 0 0 0]

would it change anything to your idea of proceeding above?

[Torsten]

Since you can set bounds on the parameters, you can also use the below code which might be more intuitive.

But the bad fit remains the same.

I'd test some initial values for the three parameters and see whether the final parameters strongly depend on the initial guess. And why do you choose 4 parameters for p0 ? After the change, you work with 3 parameters.

load frequenzee.mat
load parteImag.mat
load partReale.mat
f = frequenzee.';
epr = partReale.' %real part of your data row vector;
epi = parteImag.'; %imaginary part of your data row vector;
eps_v = [epr,epi]
%Initial values for parameters to be fitted
eps_r0 = 1;
vfw0 = 0.5;
sigma0 = 1;
%options = optimset('Display','iter','FunValCheck','on', ...
%               'MaxFunEvals',Inf,'MaxIter',Inf, ...
%               'TolFun',1e-6,'TolX',1e-6);
paramslb = [1.5 0 0];  % lower bound
paramsub = [4.5 1 Inf ];
p0 = [eps_r0,vfw0,sigma0];
% Call nonlinear optimizer
p = lsqcurvefit(@fun,p0,f,eps_v,paramslb,paramsub);%, options)
eps_r = p(1)
vfw = p(2)
vb = 1-p(2)
sigma = p(3)
stima= eps_r + vfw.*(4.9+75./(1+1i*f/18)-1i*18.*sigma./f)+...
       vb.*(2.9+55./(1+(1i*f/18).^(0.5)));
plot(f,real(stima))
hold on
plot(f,epr)
figure
plot(f,imag(stima))
hold on
plot(f,epi)
function res = fun(p,f)
  eps_r = p(1);
  vfw = p(2);
  vb = 1-p(2);
  sigma = p(3);
  eps_v_model = eps_r + vfw.*(4.9+75./(1+1i*f/18)-1i*18.*sigma./f)+...
                vb.*(2.9+55./(1+(1i*f/18).^(0.5)));  
  res = [real(eps_v_model),imag(eps_v_model)];
end  

[daniele gallo]

I work with 4 parameters because in the reference model eps-r it is assumed that it varies between [1,5-4,5]. Since we are working with new experimental data, we cannot a priori consider that this limit is respected, so it is assumed as a parameter to be estimated.

[Torsten]

eps_r is still part of the parameter set, vb is replaced by 1-vfw.