Trying to plot I-V and P-V characteristics (Am I using fsolve incorrectly?)

I am trying to plot the PV and IV characteristics for a solar array using the following code and equations. When I solve for the output current iteratively using fsolve, my current remains constant and does not "fall off" as is the expected characterisitc of a solar cell once maximum voltage is reached. Same thing with the PV curve.

I also attached the papers I used to get the equations.

%Code to generate I-V and P-V curve for one solar cell
clc
clear all 

Irradiance = 200:200:1000; %Solar Irradiance Vector in Watts/m^2
Np  = 1;       %Number of parallel panels
Ns  = 60;      %Number of series panels 
Ki  = 0.00023;  %Cell's short-circuit current temperature coefficient 
Tc  = 25+273 ;     %Cell's working temperature in Kelvin
Tref= 25+273    ;  %Cell's reference temperature in Kelvin
q   = 1.6e-19; %Electron Charge 
k   = 1.38e-23; %Boltzmann's temperature constant 
Rs  = 0.3;
n   = 2.15;     %Diode's ideal factor 
Eg  = 1.166;   %Band-gap energy 
Voc = 37.191;  %Cell's Open Circuit Voltage 
Isc = 8.449;       %Short Circuit current at 25°C and 1kW/m^2
V = linspace(0,Voc); %points from 0 to the open circuit voltage


%Calculate the reverse saturation current 
Irs = Isc./( exp((q * Voc)/(Ns * k * n * Tc)) - 1);

%Store current output Ipv and diode current Id in vectors for population
Id_array = zeros(length(Irradiance),length(V));
Iph_array = zeros(length(Irradiance),length(V));

%Create meshgrid for the voltage 
[Ir_m,V_m] = meshgrid(Irradiance,V);
Ir_m = Ir_m';
V_m = V_m';

    
for i = 1:length(V)
    %Photon Current 
    Iph = (Isc + Ki*(Tc - Tref)) * Ir_m/1000;
    Iph_array(:) = Iph;
    
    %Make current through diode a function Id
    Id =@(I) (Irs .* (exp( (q.*(V_m+I*Rs)) ./ (n * k* Tc* Ns) - 1)));
    initialguess = 0.001;
    [Id_sol] = fsolve(Id,initialguess,optimset('Display','off'));
    Id_array(:) = Id_sol;
     
     
    Ipv = Iph_array - Id_array;
     
    %Output Power
    Pout = Ipv .* V_m;
end 

    %Plot I-V and P-V curves
    subplot(2,1,1);
    plot(V,Ipv)
    xlabel('Vout'); 
    ylabel('Iout'); 

    subplot(2,1,2); 
    plot(V,Pout)
    xlabel('Vout');
    ylabel('Pout'); 




%Add legends to both subplots
subplot(2,1,1); legend('200','400','600','800','1000'); % Show legend for I-V curve
subplot(2,1,2); legend('200','400','600','800','1000'); % Show legend for P-V curve
%hold off;

[Torsten]

Your function Id you use in "fsolve" returns an array of size 5x100, but it should only return one scalar value. So something must be wrong ...

Irradiance = 200:200:1000; %Solar Irradiance Vector in Watts/m^2
Np  = 1;       %Number of parallel panels
Ns  = 60;      %Number of series panels 
Ki  = 0.00023;  %Cell's short-circuit current temperature coefficient 
Tc  = 25+273 ;     %Cell's working temperature in Kelvin
Tref= 25+273    ;  %Cell's reference temperature in Kelvin
q   = 1.6e-19; %Electron Charge 
k   = 1.38e-23; %Boltzmann's temperature constant 
Rs  = 0.3;
n   = 2.15;     %Diode's ideal factor 
Eg  = 1.166;   %Band-gap energy 
Voc = 37.191;  %Cell's Open Circuit Voltage 
Isc = 8.449;       %Short Circuit current at 25°C and 1kW/m^2
V = linspace(0,Voc); %points from 0 to the open circuit voltage


%Calculate the reverse saturation current 
Irs = Isc./( exp((q * Voc)/(Ns * k * n * Tc)) - 1);

%Store current output Ipv and diode current Id in vectors for population
Id_array = zeros(length(Irradiance),length(V));
Iph_array = zeros(length(Irradiance),length(V));

%Create meshgrid for the voltage 
[Ir_m,V_m] = meshgrid(Irradiance,V);
Ir_m = Ir_m';
V_m = V_m';
Id =@(I) (Irs .* (exp( (q.*(V_m+I*Rs)) ./ (n * k* Tc* Ns) - 1)));
initialguess = 0.001;
Id(initialguess)    

ans = 5×100

0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0004 0.0004 0.0005 0.0005 0.0006 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0004 0.0004 0.0005 0.0005 0.0006 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0004 0.0004 0.0005 0.0005 0.0006 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0004 0.0004 0.0005 0.0005 0.0006 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0004 0.0004 0.0005 0.0005 0.0006 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 1

[Arnav]

Hi Jordina Pierre,

I understand that you are trying to find the characteristic P-V and I-V plots of a PV Array. I found 2 issues in the code snippet you provided. These are shown below:

• Dimensional Mismatch:

Ir_m and V_m are 2 matrices that have dimensions 5x100. Each row corresponds to an intensity and each column corresponds to a voltage.

Since the for-loop iterates over voltages, we need to solve for a 5x1 current vector where each element corresponds to a different illumination intensity.

• Current Equation: The equation to be solved is:

In other words, we require a value of I that satisfies this equation for a particular V. fsolve function is used to find the zeros of a function. Therefore, the function handle should be this instead:

According to these points, for-loop can be modified as shown below:

for i = 1:length(V) 
    % Photon Current for all irradiance levels 
    Iph = (Isc + Ki * (Tc - Tref)) * Ir_m(:, i) / 1000; 
    Iph_array(:, i) = Iph; 
    
    % Diode current calculation for all irradiance levels 
    Id = @(I) Irs * (exp((q * (V(i) + I * Rs)) / (n * k * Tc * Ns)) - 1) - I; 
    initialguess = zeros(size(Irradiance)) + 0.001; 
    initialguess = initialguess'; 
    Id_sol = fsolve(Id, initialguess, opts); 
    Id_array(:, i) = Id_sol; 
    
    % Calculate Ipv and Pout for all irradiance levels 
    Ipv(:, i) = Iph_array(:, i) - Id_array(:, i); 
    Pout(:, i) = Ipv(:, i) * V(i); 
end 

The resultant plots are shown below:

You can learn more about fsolve function here:

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

Hope it helps!