I am working on solve switching angle of 9-level OHSW by using the Newton-Raphson method in MATLAB, actually issue is when I run for 11 level, the solution approaches to some stable value and when for 9 it does not approach to some stable solutions, the code it's work but there is small error in the coding, please any one can help me to find the error?

function p_send = nine_level_all_outputs (m_p)

close all;

E=4; %The number of dc sources
dM=0.0005; %The modulation index step
Mstart=.5; %The initial modulation index
M=Mstart*E;
Mrange=1000; %Range of calculation results

p1=20*pi/180;
p2=30*pi/180;
p3=50*pi/180;
p4=70*pi/180;
% Guess initial value of switching angles

p=[p1 p2 p3 p4]'; %The switching angle matrix

for j=1:Mrange
t=[M*pi/4 0 0 0]'; %The corresponding harmonic amplitude matrix
df=1; i=1;
    while abs(df) > 1e-20 & i < 100 %i=40 %Degree of accuracy condition
    p1=p(1,:); p2=p(2,:); p3=p(3,:); p4=p(4,:);
    f=[cos(p1)+cos(p2)+cos(p3)+cos(p4);% The nonlinear system matrix
        cos(3*p1)+cos(3*p2)+cos(3*p3)+cos(3*p4);
        cos(5*p1)+cos(5*p2)+cos(5*p3)+cos(5*p4);
        cos(7*p1)+cos(7*p2)+cos(7*p3)+cos(7*p4)];
    delf=[-sin(p1) -sin(p2) -sin(p3) -sin(p4);
        -3*sin(3*p1) -3*sin(3*p2) -3*sin(3*p3) -3*sin(3*p4);
        -5*sin(5*p1) -5*sin(5*p2) -5*sin(5*p3) -5*sin(5*p4);
        -7*sin(7*p1) -7*sin(7*p2) -7*sin(7*p3) -7*sin(7*p4)];

    df=inv(delf)*(t-f); % Calculate solution error
    p=p+df; % Update the solutions.
    p = acos(abs(cos(p)));
    i=i+1;
    end

p = sort(p);
pf1(j) = p(1);
pf2(j) = p(2);
pf3(j) = p(3);
pf4(j) = p(4);
mm(j)=Mstart; % The modulation index
Mstart = Mstart +dM;
M=Mstart*E; %Update the modulation index

end

v = floor((m_p-0.5)/0.0005)+1;
if v==1001
    v=v-1;
end

mm(v);
p_send = [pf1(v) pf2(v) pf3(v) pf4(v)]*180/pi ;
%--data for plotting--
figure(1);
pf1=abs(pf1); pf2=abs(pf2); pf3=abs(pf3); pf4=abs(pf4); %Convert neg angle to pos angle
plot(mm,pf1*180/pi,'^r',mm,pf2*180/pi,'^g',mm,pf3*180/pi,'^b',mm,pf4*180/pi,'^m');
legend('Alpha 1','Alpha 2','Alpha 3','Alpha 4')
axis([mm(1) mm(Mrange) 0 90]);

%--calculate phase THD--
for i=1:200 % Empty the harmonic matrix
ha(i)=0;
hindex(i)=i;
end

for j=1:length(pf1)
sumh=0;
sumhll=0;
    for i=1:100
        k=2*i-1; % k is odd number
        ha(k)=(cos(k*pf1(j))+cos(k*pf2(j))+cos(k*pf3(j))+cos(k*pf4(j)))/k;
% Calculate amplitude of harmonic amplitude
    end
    for i=2:200
        sumh=ha(i)^2+sumh;
    end
THDph(j)=sqrt(sumh)/ha(1)*100; % Calculate phase voltage THD


end

figure(2);
plot(mm,THDph,'k');
legend('THD of 9-level inverter');
title('The Voltage THD vs. The Modulation Index of 9-level OSHW');
xlabel('Modulation Index');
ylabel('THD (%)');
% End of File

 

this is example same the code

 

please delete this thread