Trouble using lsqnonlin to fit the IV curve of a solar cell

I am trying to fit some current density-voltage data of a solar cell which i have measured in the lab, to a non ideal solar cell equation. Basically, i have some arrays of V (voltage), and J (current density) from a solar cell measurement. I am trying to fit and extract some parameters from the equation (in this case, n, Rs and Rsh), where the equation also has some other constants and measured values which i know (namely, k, T, q, Jsc, Voc, V). This is J_cal below (basically for those interested, it is the non ideal solar cell equation rearranged, using the Lambert W function to solve for it).

I am trying to use lsqnonlin to fit the curve, to extract Rs, Rsh and n:

% % Define the initial guess for the parameters to be extracted
n0 = 2;
Rs0 = 2;
Rsh0= 500;

p0 = [n0, Rs0, Rsh0];

% Define the lower and upper bounds for the parameters to be extracted
lb = [0, 0.1, 25];
ub = [5, 100, 10000];

% Use the lsqnonlin function to find the best-fit parameters
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','Display','iter-detailed');
p = lsqnonlin( @ShockleyQueisser, p0, lb, ub, options)

and I also have the model function:

% Define the model function 
function J_model = ShockleyQueisser(p)
    global V Voc Jsc T k q J;
     n = p(1);
     Rs = p(2);
     Rsh = p(3);
    J_cal = ((n*k*T)/(q*Rs))*lambertw((q*Rs)/(n*k*T)*(Jsc-Voc/(Rs+Rsh))*(exp((-q*Voc)/(n*k*T)))*(exp(q/(n*k*T)*(Rs*Jsc+((Rsh.*V)./(Rsh+Rs))))))+V./Rs-Jsc-((Rsh.*V)./(Rs*(Rsh+Rs)));
    J_model = J_cal-J;
end

which is what is trying to be solved.

My problem is, is that it just doesn't iterate and doesn't seem to find the minimum of the sum of squares. It runs once, returns p(0) as the starting guess, and then stops:

First-Order Norm of

Iteration Func-count Residual optimality Lambda step

0 4 0 0 0.01

Optimization completed: The final point is the initial point.

The first-order optimality measure, 0.000000e+00, is less than

1e-4*options.FunctionTolerance = 1.000000e-10.

I am confident J_cal is correct, since if i paste the equation into the command window after i set the initial guess, it gives me something which looks like a good starting guess, but my problem is the code wont iterate.

The global variables are all 1x1, apart from V and J which are 91x1.

I think the issue is how p is defined... in the online documentation, it says i just define the function (in this case 'ShockleyQueisser') however for some reason other help points at setting a handle, @, and i don't really know why. Setting the handle this way makes the code run, but iterates only once. Not putting it gives an error.

[Chris]

Could you please provide some representative values for the globals?

[Chris]

The exit condition is indicating J==J_cal with your initial guess and, since your objective function is the difference between the two, the problem is minimized. You can't get any smaller than zero.

Really though, you could go negative, which is why we minimize the sum of squares--in this case, you would want:

J_model = (J_cal-J).^2;

Here I'm assuming J is what you measured and J_cal is your theoretical model. Because you're handing the solver J_model as your objective function to be minimized.

[Torsten]

@Chris

For lsqnonlin,

sum_i f(i)^2

is minimized, and the user is asked to provide the f(i).

Thus returning

res = J_calc - J_experiment

is correct (if J_experiment is an array of experimental measurements of the same size as J_calc).

Is it possible that the values returned by lambertW are complex and that this is the reason lsqnonlin stops ?

[Jake Bowers]

Hello both. @Chris, below is the full script running, with an attachment of a data text file which is the data i am trying to fit (plus also tidy up_. It provides the values of Jsc, Voc, plus the measured data J and V.

The first guess it does is not near the final solution. You can change it and see and compare against the measured data.

@Torsten, i dont believe they are complex. As mentioned above, if i copy the equation in J_cal into the command window, and run it, i get a series of OK (for a first guess) data. No complex numbers there.

clc 
clear

%reads input file and extracts the IV parameters taken from Labview
%generated txt file
inputfile='data.txt'; %sets input file name

IV_param=readtable(inputfile, 'Headerlines',3, 'ReadRowNames', true);
IV_param(21:end,:)=[];
IV_param(:,1)=[];

%reads input file and extracts I and V data and puts into an array from a
%table
IV = readtable(inputfile, 'HeaderLines', 24);
IV=table2array(IV);

%Trims data to remove anything above the compliance limit 
IV(IV(:,2)>= abs(IV_param.Var2(8,:)),:)=[]; 

% Sets voltage and current arrays
V=IV(:,1);
I=IV(:,2);

%Calculates current density in mA/cm2 from Labview defined area
J=(IV(:,2)/IV_param.Var2(4,:));

%constants 
k = 1.380649e-23; % Boltzmann constant
T = 293; %Temp in Kelvin
q = 1.6e-19; %charge of an electron.

%Interpolate to estimate Voc
Voc=interp1(J,V,0);

%Interpolate to estimate Jsc
Jsc=interp1(V,J,0);
Jsc=-Jsc;

% % Define the initial guess for the parameters to be extracted
n0 = 2;
Rs0 = 2;
Rsh0= 500;

p0 = [n0, Rs0, Rsh0];

% Define the lower and upper bounds for the parameters to be extracted
lb = [0, 0.1, 25];
ub = [5, 100, 10000];

% Use the lsqnonlin function to find the best-fit parameters
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','Display','iter-detailed','MaxIterations',5000);
[p,exitflag] = lsqnonlin( @(p) ShockleyQueisser(p), p0, lb, ub, options);

%Defines Rs, Rsh and n in p
Rs = p(2);
Rsh = p(3);
n = p(1);

%Plot the original J-V data
figure;
plot(V, J, 'o');
xlabel('Voltage (V)');
ylabel('Current Density (mA/cm^2)');
legend('Experimental Data');

% Define the model function 
function J_model = ShockleyQueisser(p)
    global V Voc Jsc T k q J;
     n = p(1);
     Rs = p(2);
     Rsh = p(3);
    J_cal = ((n*k*T)/(q*Rs))*lambertw((q*Rs)/(n*k*T)*(Jsc-Voc/(Rs+Rsh))*(exp((-q*Voc)/(n*k*T)))*(exp(q/(n*k*T)*(Rs*Jsc+((Rsh*V)/(Rsh+Rs))))))+V/Rs-Jsc-((Rsh*V)/(Rs*(Rsh+Rs)));
     %J_cal = ((n*k*T)/(q*Rs))*lambertw((q*Rs)/(n*k*T)*(Jsc-Voc/(Rs+Rsh))*(exp((-q*Voc)/(n*k*T)))*(exp(q/(n*k*T)*(Rs*Jsc+((Rsh.*V)./(Rsh+Rs))))))+V./Rs-Jsc-((Rsh.*V)./(Rs*(Rsh+Rs)));
    J_model = J_cal-J;
end

Answer 1

[Torsten]

You did not include

global V Voc Jsc T k q J

in the script part.

clc 
clear

global V Voc Jsc T k q J
%reads input file and extracts the IV parameters taken from Labview
%generated txt file
inputfile='data.txt'; %sets input file name

IV_param=readtable(inputfile, 'Headerlines',3, 'ReadRowNames', true);
IV_param(21:end,:)=[];
IV_param(:,1)=[];

%reads input file and extracts I and V data and puts into an array from a
%table
IV = readtable(inputfile, 'HeaderLines', 24);
IV=table2array(IV);

%Trims data to remove anything above the compliance limit 
IV(IV(:,2)>= abs(IV_param.Var2(8,:)),:)=[]; 

% Sets voltage and current arrays
V=IV(:,1);
I=IV(:,2);

%Calculates current density in mA/cm2 from Labview defined area
J=(IV(:,2)/IV_param.Var2(4,:));

%constants 
k = 1.380649e-23; % Boltzmann constant
T = 293; %Temp in Kelvin
q = 1.6e-19; %charge of an electron.

%Interpolate to estimate Voc
Voc=interp1(J,V,0);

%Interpolate to estimate Jsc
Jsc=interp1(V,J,0);
Jsc=-Jsc;

% % Define the initial guess for the parameters to be extracted
n0 = 2;
Rs0 = 2;
Rsh0= 500;

p0 = [n0, Rs0, Rsh0]

p0 = 1×3

2 2 500

% Define the lower and upper bounds for the parameters to be extracted
lb = [0, 0.1, 25];
ub = [5, 100, 10000];

% Use the lsqnonlin function to find the best-fit parameters
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','Display','iter-detailed','MaxIterations',5000);
[p,exitflag] = lsqnonlin( @(p) ShockleyQueisser(p), p0, lb, ub, options);

                                        First-Order                    Norm of 
 Iteration  Func-count    Residual       optimality      Lambda           step
     0           4      0.00147924         0.00116         0.01
     1           8      0.00114861         0.00113        0.001       0.121731
     2          12     8.10179e-05        0.000472       0.0001       0.512376
     3          16      2.6781e-06        2.98e-05        1e-05       0.204876
     4          20     5.02191e-07        1.53e-06        1e-06       0.126386
     5          24      4.8159e-07        3.37e-08        1e-07      0.0143882
     6          28     4.81585e-07        3.22e-10        1e-08     0.00252563
     7          32     4.81572e-07        2.51e-10        1e-09      0.0251263
     8          36     4.81451e-07        2.44e-10        1e-10       0.244337
     9          40     4.80626e-07        3.46e-10        1e-11        1.90361
    10          44     4.79421e-07        5.89e-09        1e-12        5.01446
    11          48     4.79319e-07        4.77e-10        1e-13        1.87756
    12          52     4.79319e-07        8.71e-11        1e-14      0.0981056

Optimization completed: The relative first-order optimality measure, 8.707231e-11,
is less than 1e-4*options.FunctionTolerance = 1.000000e-10.

J_model = ShockleyQueisser(p);
J_cal = J_model + J;
hold on
plot(V,J_cal)
%Defines Rs, Rsh and n in p
n = p(1)

n = 1.8934

Rs = p(2)

Rs = 1.1913

Rsh = p(3)

Rsh = 509.1659

%Plot the original J-V data
%figure;
plot(V, J, 'o');
xlabel('Voltage (V)');
ylabel('Current Density (mA/cm^2)');
legend('Experimental Data');
hold off

% Define the model function 
function J_model = ShockleyQueisser(p)
    global V Voc Jsc T k q J
     n = p(1);
     Rs = p(2);
     Rsh = p(3);
    J_cal = ((n*k*T)/(q*Rs))*lambertw((q*Rs)/(n*k*T)*(Jsc-Voc/(Rs+Rsh))*(exp((-q*Voc)/(n*k*T)))*(exp(q/(n*k*T)*(Rs*Jsc+((Rsh*V)/(Rsh+Rs))))))+V/Rs-Jsc-((Rsh*V)/(Rs*(Rsh+Rs)));
     %J_cal = ((n*k*T)/(q*Rs))*lambertw((q*Rs)/(n*k*T)*(Jsc-Voc/(Rs+Rsh))*(exp((-q*Voc)/(n*k*T)))*(exp(q/(n*k*T)*(Rs*Jsc+((Rsh.*V)./(Rsh+Rs))))))+V./Rs-Jsc-((Rsh.*V)./(Rs*(Rsh+Rs)));
    J_model = J_cal-J;
end

[Jake Bowers]

Wow! Do I feel stupid... did the hard bit, and forgot the easy thing at the beginning!

Thank you so much for your help... you made my week!