Problem in getting correct output for double diode model of pv array

the following is the code for solving current equation of a double diode model of a pv array, it throws an error when I run it

% Parameters

I_L = 8.3; % Light-generated current

I_01 = 1.22e-12; % Reverse saturation current for diode 1

I_02 = 1.22e-12; % Reverse saturation current for diode 2

n_1 = 1.36; % Ideality factor for diode 1

n_2 = 2.34; % Ideality factor for diode 2

k = 1.380649e-23; % Boltzmann constant

q = 1.60217662e-19; % Electronic charge

T = 298.15; % Temperature

R_S = 0.221; % Series resistance

% Voltage range

V_min = 0;

V_max = 32.9;

V_step = 0.01;

V = V_min:V_step:V_max;

% Double diode model

I = zeros(size(V));

for i = 1:length(V)

f = @(x) I_L - I_01*(exp(q*(V(i)+x*R_S)/(n_1*k*T))-1) - I_02*(exp(q*(V(i)+x*R_S)/(n_2*k*T))-1) - x/R_S;

I(i) = fsolve(f, 33);

end

% Plot current vs voltage

figure;

plot(V,I);

xlabel('Voltage (V)');

ylabel('Current (A)');

title('Double Diode Model');

[Aryan Sharma]

The correct output should look like this

[Vinayak]

When I run the following code snippet on my machine,

I receive an error pertaining to the fsolve function. The error suggests that the function being passed to the fsolve solver is returning undefined values for the initial point used as an argument in the solver. This is likely caused by a discrepancy in one of the input parameters or an incorrect set of initial values.

[Alan Stevens]

The last term in your function (x/R_S) seems to have units of current/resistance (or would do if units were explicitly included!), whereas the other three terms have units of current only.

[Joel Van Sickel]

It's likely your equation for fsolve is incorrect, can you share the analytical formula that you are tryint to solve?

[Aryan Sharma]

The equations used in the code are as follows:

• Equation (1): Light-generated current (IL):
• IL = I_L - I_01 * (exp(q * (V(i) + x * R_S) / (n_1 * k * T)) - 1) - I_02 * (exp(q * (V(i) + x * R_S) / (n_2 * k * T)) - 1)

• Equation (2): Diode current for diode 1 (ID1):
• ID1 = I_01 * (exp(q * (V(i) + x * R_S) / (n_1 * k * T)) - 1)

• Equation (3): Diode current for diode 2 (ID2):
• ID2 = I_02 * (exp(q * (V(i) + x * R_S) / (n_2 * k * T)) - 1)

• Equation (4): Equation to solve for the current (I) at each voltage point:
• f = @(x) I_L - I_01 * (exp(q * (V(i) + x * R_S) / (n_1 * k * T)) - 1) - I_02 * (exp(q * (V(i) + x * R_S) / (n_2 * k * T)) - 1) - x / R_S

These equations model the behavior of the double-diode photovoltaic module and are used to calculate the current (I) at different voltage (V) points within the specified range. The fsolve function is utilized to numerically solve the equation (4) and obtain the corresponding current values.

Answer 1

[Sayan]

Hi @Aryan Sharma,

I understand from your issue that the code you are trying to run to obtain the "I-V" characteristics of a PV cell array is throwing an error. Below are the possible workarounds for the issue:

• The number of function evaluations can be increased for the "fsolve" function to overcome the error.
• The initial condition to determine the value of current "I" defined in the "fsolve" function can be provided to the previous computed value of current rather than hard coding it to 33.
• The last term in the function handle expression "f" should be dictating the current through the shunt resistance. But shunt resistance is not declared in the code, and the expression "x/R_s" is wrong.
• The changes are shown in the below-mentioned code snippet. The code mentioned below should be included after the line "I = zeros(size(V));"

%declare the value of shunt resistance 
R_sh=100;%say 100
%compute the initial value of current for 0 voltage
f1 = @(x) I_L - I_01*(exp(q*(x*R_S)/(n_1*k*T))-1) - I_02*(exp(q*(x*R_S)/(n_2*k*T))-1) - x*R_S/R_sh;
I(1)=fsolve(f1,0);
%declare the number of function evaluations for "fsolve" function
options = optimoptions('fsolve', 'MaxFunctionEvaluations', 1000);
for i = 2:length(V)
    f = @(x) I_L - I_01*(exp(q*(V(i)+x*R_S)/(n_1*k*T))-1) - I_02*(exp(q*(V(i)+x*R_S)/(n_2*k*T))-1) - x*R_S/R_sh;%declare x*R_S/R_sh as the current through shunt branch
    I(i) = fsolve(f,I(i-1),options);%provide the "options" to "fsolve" function call
end
%now plot the I-V data and change the value of the constants to meet the
%required curve 

Further information on "fsolve" function can be found in the following documentation.

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

Hope this helps in resolving the issue.