Determining the series and parallel resistances of a two diode PV model in matlab by iteration

Hello peeps ! I am currently on an internship and my topic requires me to simulate PV cells in MATLAB. I am using a two diode model, and I have written a code to determine Rs and Rp by iteration, but the problem is the code kept running for more than 14 hours ( I had to stop it, it did not stop by itself) and I did not obtain any results. Here is the code:

clear;
clc;
q=1.6*10^(-19);%electron charge
K=1.38064852*10^(-23);%Boltzmann constant
NOCT=45;%the nominal operating cell temperature
Vmpp=30.12;%the voltage at the maximum power point 
Impp=8.3;%the current at the maxium power point
Voc= 37.85;%the open circuit voltage
Isc=8.65;%the short circuit current
Ki=0.005;%temperature coefficient of Isc
Kv=-0.17;% temperature coefficient of Voc
Nc=60;% number of cells connected in series in the panel
Tmean=26.45;% Temperature Annual Daily Mean
Emean=425.78;% Radiation Annual Daily Mean
Pmaxe=Vmpp*Impp;% experimental maximum power point
Rs=0;%initializing the series resistance
Tjpn=273.15+Tmean+(NOCT-20)*(Emean/800);% calculating the junction temperature
Ipv=Isc+Ki*(Tjpn-298.15)*(Emean/1000);%the current generated by the incidence of light
Isat=(Isc+Ki*(Tjpn-298.15))/(exp((q*(Voc+(Kv*(Tjpn-298.15))))/(Nc*K*Tjpn))-1);%the diode saturation current
Rpmin= Vmpp/(Ipv-2*Isat*exp((q*Vmpp)/(Nc*K*Tjpn)-1)-Impp);%initializing the parallel resistance
Rp=Rpmin;
error=1;
while error>0.05

    Rp=(Vmpp+Impp*Rs)/(Ipv-2*Isat*exp((q*(Vmpp+Impp*Rs)/(Nc*K*Tjpn))-1)-Impp);
    %calculating teh I-V equation
    clear V;
    clear I;
    V=0:0.01:38;
    I=zeros(1,size(V,2));
    g=zeros(1,size(V,2));
    d=zeros(1,size(V,2));
    I_=zeros(1,size(V,2));
    for j=1:size(V)% calculating for all voltage values
        g(j)=Ipv-2*Isat*exp((q*(V(j)+I(j)*Rs))/(Nc*K*Tjpn))-(V(j)+I(j)*Rs)/Rp-I(j);
        %solving g=I-f(I,V)=0 using Newton-Raphson method
        while(abs(g(j))>0.05)
        g(j)=Ipv-2*Isat*exp((q*(V(j)+I(j)*Rs))/(Nc*K*Tjpn))-(V(j)+I(j)*Rs)/Rp-I(j);
        d(j)=2*Isat*((q*Rs)/(Nc*K*Tjpn))*exp((q*(V(j)+I(j)*Rs))/(Nc*K*Tjpn))-(Rs/Rp)-1;%calculating the derivative of g(j)
        I_(j)=I(j)-g(j)/d(j);
        I(j)=I_(j);

        end %while(abs(g(j))>tol) 
    end %for j=1:size(V,2)
        P=(Ipv-2*Isat*exp((q*(V+I.*Rs))/(Nc*K*Tjpn))-(V+I.*Rs)/Rp).*V;
        Pmaxm=max(P);
        error=abs(Pmaxe-Pmaxm);%checking if the calculated Pmax is closer to the tabulated one on the manufacturer datasheet
        Rs=Rs+0.01;%incrementing Rs
end %while error>tol
fprintf('Rs is %g and Rp is %g',Rs,Rp)
In order to evaluate the output current I used Newton-Raphson method.

Can you please help me finding the error in my code... I really need your help guys to proceed with my internshship's topic.

Thank you in advance :)

[David Goodmanson]

Hi Wissal, I am not commenting on the code overall and only checked out the NR part. There appears to be a sign error between g(j) and d(j) in the 2*Isat*exp(....) term which can't be helping.