Fitting routine, optimization with two functions for one problem

Hello matlab team. I have a question to an optimization of a fitting routine. I am trying to plot data to a set of functions. The problem is that I do not know how to handle two functions for fitting one problem. At the moment I am trying to do the following with the lsqcurvefit. Data attached.

clear all;
close all;

%%Prepare data
% Load data
Aktor = importdata(['Aktor1.txt']);
Aktor_freq = Aktor(:,1);
Aktor_ampl = Aktor(:,2);
Aktor_phase = Aktor(:,3);

Re_Y = Aktor_ampl.*cos(deg2rad(Aktor_phase));   % Conductance G
Im_Y = Aktor_ampl.*sin(deg2rad(Aktor_phase));   % Susceptance B

plot(Re_Y,Im_Y,'ok','LineWidth',1.5);            % Show data in complex plain
hold all

% Define start parameters
C_0 = 3.8182e-08;
R_m = 8.0856e+03;
L_m = 15.3030;
C_m = 1.0271e-09;
tan_d = 0.02;

x_data =  Aktor_freq;                           % x-values
y_data =  Aktor_ampl;                           % y-values

% Show start parameters
y_0 = (2*pi*x_data*C_0*(i+tan_d))+1./(R_m+i*(2*pi*x_data*L_m-1./(2*pi*x_data*C_m)));
plot(real(y_0),imag(y_0),'r','LineWidth',1.5);

%%Generate Fit 
%    [C_0,  R_m,  L_m,  C_m,  tan_d];
x0 = [1,    1,    1,    1,    1];                     % Relative altering 
lb = [1/10, 1/10, 1/10, 1/10, 1/10];                  % Lower bound
ub = [1*10, 1*10, 1*10, 1*10, 1*10];                  % Upper bound

Fun = @(x,xdata) abs((2*pi*xdata*C_0*x(1)*(i+tan_d*x(5)))+1./(x(2)*R_m+i*(2*pi*xdata*L_m*x(3)-1./(2*pi*xdata*C_m*x(4)))));
options = optimset('TolX',1e-20,'TolFun',1e-20);
x = lsqcurvefit(Fun,x0,x_data,y_data,lb,ub,options);

% Show fitted data
 y_fit1 = (2*pi*x_data*C_0*x(1)*(i+tan_d*x(5)))+1./(x(2)*R_m+i*(2*pi*x_data*L_m*x(3)-1./(2*pi*x_data*C_m*x(4))));
 plot(real(y_fit1),imag(y_fit1),'b','LineWidth',1.5);
  xlabel('real part')
  ylabel('imaginary part')

You can see. I am using only one function with five parameters that I fit relativ to its start value. The fit function is not that good. Showing the data in the complex plain is required. The curves does not match properly. I get a very nice result if I change tan_d = 0.05;

I would like to change the routine. I have two different functions.

w = 2*pi*Aktor_freq;

The real part is defined by:

Fun1 = R_m./(R_m^2 + (w*L_m - 1./(w*C_m)).^2) + w*tan_d*C_0 == Re_Y

The imaginary part is defined by:

Fun2 = (-w*L_m + 1./(w*C_m))./(R_m^2 + (w*L_m - 1./(w*C_m)).^2) + w*C_0 == Im_Y

So I show you the figures of the imaginary and the real part.

Real

tan_d = 0.2

Imaginary

The imaginary part seems to fit quite nice but the real part doesn't hit the spot. I found out that the global peak of the real part and the peak of Fun1 has to be almost equal.

Maybe I can somehow try to give the points around the peak of the real part a higher weight. So the fit converges to it more then to the other points. The main idea is to find the parameters so that the fit matches the ellipse of the data in the complex plain.

Is there any way to solve this problem?. Imagine how to handle both functions for a fitting procedure is out of my level. I appreciate any advice. I tried working with the genetic algorithm but with no succes.

Regards Henrik

[Henrik Schädlich]

At the moment I am trying to do the following. I want to solve the problem with fsolve documentation. It shoud be somehow possible.

First I define my functions.

function F = nonlinmodel( x, C_0, R_m, L_m, C_m, tan_d, w, Re_Y, Im_Y)
%       C_0, R_m, L_m, C_m, tan_d are the parameters I have to fit.
%
%       x(1) to x(5) are there relativ changes
%
%       w is the step of the x-Axis
%
%       Re_Y and Im_Y are the data sets.
%
%       F(1) and F(2) are almost an array of 0. As in the documentation
%           They are not exactly 0 due to the error
%           in the inital condition

% Real part
% Function out of theory minus measured data
F(1) = x(2)*R_m./((x(2)*R_m).^2 + (w*x(3)*L_m - 1./(w*x(4)*C_m)).^2) + w*x(5)*tan_d*x(1)*C_0 - Re_Y;      

% Imaginary part
% Function out of theory minus measured data
F(2) = (-w*x(3)*L_m + 1./(w*x(4)*C_m))./((x(1)*R_m).^2 + (w*x(3)*L_m - 1./(w*x(4)*C_m)).^2) + w*x(1)*C_0 - Im_Y;
end

Afterwards I optimize my inital conditions. At least I try. ^^

clear all
close all

%%Load data
Aktor = importdata(['Aktor1.w_l']);
Aktor_freq = Aktor(:,1);
Aktor_ampl = Aktor(:,2);
Aktor_phase = Aktor(:,3);

Re_Y = Aktor_ampl.*cos(deg2rad(Aktor_phase));
Im_Y = Aktor_ampl.*sin(deg2rad(Aktor_phase));
w = 2*pi*Aktor_freq;

%%Define inital values
C_0 = 3.948675967565607e-08;
R_m = 7.120051155668139e+03;
L_m = 25.35529039370214e+00;
C_m = 6.242970891780357e-10;
tan_d = 0.054859859155354e-00;

%%relativ change of the inital values
x1 = [1,1,1,1,1];

%%Start fit routine
f = @(x)nonlinmodel(x, C_0, R_m, L_m, C_m, tan_d, w, Re_Y, Im_Y);
optimal_x = fsolve(f,x1);

After that I am getting an error message.

Subscripted assignment dimension mismatch.

Error in nonlinmodel (line 16)
F(1) = x(2)*R_m./((x(2)*R_m).^2 + (w*x(3)*L_m - 1./(w*x(4)*C_m)).^2) +
w*x(5)*tan_d*x(1)*C_0 - Re_Y;

Error in Fitting>@(x)nonlinmodel(x,C_0,R_m,L_m,C_m,tan_d,w,Re_Y,Im_Y)

Error in fsolve (line 230)
            fuser = feval(funfcn{3},x,varargin{:});

Error in Fitting (line 26)
optimal_x = fsolve(f,x1);

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

What am I doing wrong? I cannot figure that out. Can you help me? Best regards Henrik

[Alan Weiss]

The w argument is a vector of 600 entries. But F(1) should be a scalar.

If you want F to have 600 or 1200 entries, then you should write the nonlinmodel function differently, maybe the first line should be something like

F(1:600) = x(2)...

Alan Weiss

MATLAB mathematical toolbox documentation

[Henrik Schädlich]

Ok you are right. I changed the input with the figure handle. My result is:

clear all
close all

%%Load data
Aktor = importdata(['Aktor1.txt']);
Aktor_freq = Aktor(:,1);
Aktor_ampl = Aktor(:,2);
Aktor_phase = Aktor(:,3);

Re_Y = Aktor_ampl.*cos(deg2rad(Aktor_phase));
Im_Y = Aktor_ampl.*sin(deg2rad(Aktor_phase));
f = (1:1:3000)';

%%Start values
C_0 = 3.948675967565607e-08;
R_m = 7.120051155668139e+03;
L_m = 25.35529039370214e+00;
C_m = 6.242970891780357e-10;
tan_d = 0.054859859155354e-00;

%%relative change of start values
%x_0=[C_0, R_m, L_m, C_m, tan_d]
x =  [  1,   1,   1,   1,    1];
lb = [0.1, 0.1, 0.1, 0.1, 0.1,];
ub = [ 10,  10,  10,  10,   10];

%%define functions
Fun = {@(f)(x(2)*R_m./((x(2)*R_m).^2 + (2*pi*f*x(3)*L_m - 1./(2*pi*f*x(4)*C_m)).^2) + 2*pi*f*x(5)*tan_d*x(1)*C_0); ...
       @(f)(-2*pi*f*x(3)*L_m + 1./(2*pi*f*x(4)*C_m))./((x(1)*R_m).^2 + (2*pi*f*x(3)*L_m - 1./(2*pi*f*x(4)*C_m)).^2) + 2*pi*f*x(1)*C_0}; 

%%Fit routine
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective', ...
          'StepTolerance',1.0000e-10, 'FunctionTolerance', 1.0000e-10);
optimal_x = lsqnonlin (Fun,x,lb,ub,options);

%%Show results
figure(1);
set(gcf, 'Units', 'normalized', 'Position', [0, 0.5, 0.3, 0.3]);
plot(Aktor_freq, Re_Y,'*b','MarkerSize',1); hold all;
plot(f, Fun{1}(f),'r');
ylabel('Realteil');
xlabel('Frequenz');

figure(2);
set(gcf, 'Units', 'normalized', 'Position', [0.35, 0.5, 0.3, 0.3]);
plot(Aktor_freq, Im_Y, '*b','MarkerSize',1); hold all;
plot(f, Fun{2}(f),'r');
ylabel('Imaginärteil');
xlabel('Frequenz');

figure(3);
set(gcf, 'Units', 'normalized', 'Position', [0.7, 0.5, 0.3, 0.3]);
plot(Re_Y, Im_Y, 'b*','MarkerSize',1); hold all; 
plot(Fun{1}(f), Fun{2}(f),'r');
ylabel('Imaginärteil');
xlabel('Realteil');

The problem is, that the algorithm does not change the fit parameters. What do I do wrong here? Do you have any suggestions?

Regards Henrik

[Alan Weiss]

I am not sure, but you seem to be summing squared differences within your objective function. If you are, then don't do that, just pass in the vector of values to lsqnonlin. The solver does the summing and squaring internally.

Also, is your objective function a cell array of function handles or a plain function handle? I believe that it should be a plain function handle. For multidimensional output, just have your objective function return a vector or matrix, but don't have a cell array as your objective function.

You also calculate things using functions like cos(deg2rad(...)). MATLAB has functions that calculate in degrees directly, such as cosd.

Alan Weiss

MATLAB mathematical toolbox documentation