warning: badly scaled matrix

hello i have a problem about this code:

---------------------------------------- this is the main code:

-------------

    %The following convention is used for the four types of buses available
    %in conventional power flow studies:

    %bustype = 1 is slack or swing bus

    %bustype = 2 is generator PV bus

    %bustype = 3 is load PQ bus

    %bustype = 4 is generator PQ bus

    %
    %The five buses in the network shown in Figure 4.6 are numbered for the
    % purpose of the power flow solution, as follows:
    %North = 1
    %South = 2
    %Lake = 3
    %Main = 4
    %Elm = 5
    %
    %Bus data
    %nbb = number of buses
    %bustype = type of bus
    %VM = nodal voltage magnitude
    %VA = nodal voltage phase angle
    nbb = 6;
    bustype(1) = 1 ; VM(1) = 1.06 ; VA(1) =0 ;
    bustype(2) = 2 ; VM(2) = 1 ; VA(2) =0 ;
    bustype(3) = 3 ; VM(3) = 1 ; VA(3) =0 ;
    bustype(4) = 3 ; VM(4) = 1 ; VA(4) =0 ;
    bustype(5) = 3 ; VM(5) = 1 ; VA(5) =0 ;
    bustype(6) = 5 ; VM(6) = 1 ; VA(6) =0 ;
    %
    %Generator data
    %ngn = number of generators
    %genbus = generator bus number
    %PGEN = scheduled active power contributed by the generator
    %QGEN = scheduled reactive power contributed by the generator
    %QMAX = generator reactive power upper limit
    %QMIN = generator reactive power lower limit
    ngn = 2 ;
    genbus(1) = 1 ; PGEN(1) = 0 ; QGEN(1) = 0 ; QMAX(1) = 5 ; QMIN(1) = -5 ;
    genbus(2) = 2 ; PGEN(2) = 0.4 ; QGEN(2) = 0 ; QMAX(2) = 3 ; QMIN(2) = -3 ;
    %
    %Transmission line data
    %ntl = number of transmission lines
    %tlsend = sending end of transmission line
    %tlrec = receiving end of transmission line
    %tlresis = series resistance of transmission line
    %tlreac = series reactance of transmission line
    %tlcond = shunt conductance of transmission line
    %tlsuscep = shunt susceptance of transmission line
    ntl = 7 ;
    tlsend(1) = 1 ; tlrec(1) = 2 ; tlresis(1) = 0.02 ; tlreac(1) = 0.06 ;
    tlcond(1) = 0 ; tlsuscep(1) = 0.06 ;
    tlsend(2) = 1 ; tlrec(2) = 3 ; tlresis(2) = 0.08 ; tlreac(2) = 0.24 ;
    tlcond(2) = 0 ; tlsuscep(2) = 0.05 ;
    tlsend(3) = 2 ; tlrec(3) = 3 ; tlresis(3) = 0.06 ; tlreac(3) = 0.18 ;
    tlcond(3) = 0 ; tlsuscep(3) = 0.04 ;
    tlsend(4) = 2 ; tlrec(4) = 4 ; tlresis(4) = 0.06 ; tlreac(4) = 0.18 ;
    tlcond(4) = 0 ; tlsuscep(4) = 0.04 ;
    tlsend(5) = 2 ; tlrec(5) = 5 ; tlresis(5) = 0.04 ; tlreac(5) = 0.12 ;
    tlcond(5) = 0 ; tlsuscep(5) = 0.03 ;
    tlsend(6) = 6 ; tlrec(6) = 4 ; tlresis(6) = 0.01 ; tlreac(6) = 0.03 ;
    tlcond(6) = 0 ; tlsuscep(6) = 0.02 ;
    tlsend(7) = 4 ; tlrec(7) = 5 ; tlresis(7) = 0.08 ; tlreac(7) = 0.24 ;
    tlcond(7) = 0 ; tlsuscep(7) = 0.05 ;
    %
    %Shunt data
    %nsh = number of shunt elements
    %shbus = shunt element bus number
    %shresis = resistance of shunt element
    %shreac = reactance of shunt element:
    %+ve for inductive reactance and –ve for capacitive reactance
    nsh = 0 ;
    shbus(1) = 0 ; shresis(1) = 0 ; shreac(1) = 0 ;
    %
    %Load data
    %nld = number of load elements
    %loadbus = load element bus number
    %PLOAD = scheduled active power consumed at the bus
    %QLOAD = scheduled reactive power consumed at the bus
    nld = 4 ;
    loadbus(1) = 2 ; PLOAD(1) = 0.2 ; QLOAD(1) = 0.1 ;
    loadbus(2) = 3 ; PLOAD(2) = 0.45 ; QLOAD(2) = 0.15 ;
    loadbus(3) = 4 ; PLOAD(3) = 0.4 ; QLOAD(3) = 0.05 ;
    loadbus(4) = 5 ; PLOAD(4) = 0.6 ; QLOAD(4) = 0.1 ;
    itmax = 100;
    tol = 1e-12;
    nmax = 2*nbb;
    NLTC = 1 ;
    LTCsend(1) = 3 ; LTCrec(1) = 6 ; Rltc(1) = 0 ; Xltc(1) = 0.1 ;
    Tap(1) = 1 ; TapHi(1) = 1.5 ; TapLo(1) = 0.5 ; Bus(1) = 3 ; LTCVM(1) = 1 ;
    [YR,YI] = YBus(tlsend,tlrec,tlresis,tlreac,tlsuscep,tlcond,shbus,...
    shresis,shreac,ntl,nbb,nsh);
    [VM,VA,it,Tap] = LTCNewtonRaphson(tol,itmax,ngn,nld,nbb,bustype,...
    genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA,NLTC,...
    LTCsend,LTCrec,Rltc,Xltc,Tap,TapHi,TapLo,Bus,LTCVM,nmax);
    [PQsend,PQrec,PQloss,PQbus] = PQflows(nbb,ngn,ntl,nld,genbus,...
    loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,QLOAD,...
    VM,VA);
    [LTCPQsend,LTCPQrec] = LTCPQflows(NLTC,LTCsend,LTCrec,Rltc,Xltc,...
    Tap,VM,VA);
    it %Iteration number
    VM %Nodal voltage magnitude (p.u.)
    VA = VA*180/pi %Nodal voltage phase angle(Deg)
    PQsend %Sending active and reactive powers (p.u.)
    PQrec %Receiving active and reactive powers (p.u.)
    Tap %Final transformer tap position
    % End of Main LTCNewtonRaphson PROGRAM

    --------------------------

    function [VM,VA,it,Tap] = LTCNewtonRaphson(tol,itmax,ngn,nld,nbb,...
    bustype,genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,...
    VA,NLTC,LTCsend,LTCrec,Rltc,Xltc,Tap,TapHi,TapLo,Bus,LTCVM,nmax)
    % GENERAL SETTINGS
    D=zeros(1,nmax);
    flag = 0;
    it = 1;
    % CALCULATE NET POWERS
    [PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,QGEN,PLOAD,...
    QLOAD);
    while (it <= itmax && flag==0)
    % CALCULATED POWERS
    [PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI);
    [QNET,bustype] = GeneratorsLimits(ngn,genbus,bustype,QGEN,QMAX,...
    QMIN,QCAL,QNET, QLOAD, it, VM, nld, loadbus);
    % CALCULATED LTC POWERS
    [PCAL,QCAL,ltcPCAL,ltcQCAL] = LTCCalculatedPowers(NLTC,LTCsend,...
    LTCrec,Tap,Rltc,Xltc,VM,VA,PCAL,QCAL);
    % POWER MISMATCHES
    [DPQ,DP,DQ,flag] = PowerMismatches(nmax,nbb,tol,bustype,...
    flag,PNET,QNET,PCAL,QCAL);
    % Check for convergence
    if flag == 1
    break
    end
    % JACOBIAN FORMATION
    [JAC] = NewtonRaphsonJacobian(nmax,nbb,bustype,PCAL,QCAL,...
    VM,VA,YR,YI);
    % LTC JACOBIAN UPDATING
    [JAC] = LTCNewtonRaphsonJacobian(bustype,LTCsend,LTCrec,NLTC,Tap,...
    Bus,Rltc,Xltc,ltcPCAL,ltcQCAL,VM,VA,JAC);
    % SOLVE FOR THE STATE VARIABLES VECTOR
    D = JAC\DPQ';
    % UPDATE STATE VARIABLES
    [VA,VM] = StateVariablesUpdates(nbb,D,VA,VM);
    % UPDATE LTC TAPs
    [VM,Tap] = LTCUpdates(VM,D,bustype,NLTC,LTCsend,LTCrec,Tap,Bus,...
    LTCVM);
    % CHECK FOR POSSIBLE LTC TAPs’ LIMITS VIOLATIONS
    [Tap,bustype] = LTCLimits(bustype,NLTC,Tap,TapHi,TapLo,LTCsend,...
    LTCrec);
    it = it + 1;
    end

------------------------

%Function to calculate the net scheduled powers
function [PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,...
QGEN, PLOAD,QLOAD)
% CALCULATE NET POWERS
PNET = zeros(1,nbb);
QNET = zeros(1,nbb);
for ii = 1: ngn
PNET(genbus(ii)) = PNET(genbus(ii)) + PGEN(ii);
QNET(genbus(ii)) = QNET(genbus(ii)) + QGEN(ii);
end
for ii = 1: nld
PNET(loadbus(ii)) = PNET(loadbus(ii)) - PLOAD(ii);
QNET(loadbus(ii)) = QNET(loadbus(ii)) - QLOAD(ii);
end
%End function NetPowers

-------------------------------

%Function to calculate injected bus powers
function [PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI)
% Include all entries
PCAL = zeros(1,nbb);
QCAL = zeros(1,nbb);
for ii = 1: nbb
PSUM = 0;
QSUM = 0;
for jj = 1: nbb
PSUM = PSUM + VM(ii)*VM(jj)*(YR(ii,jj)*cos(VA(ii)-VA(jj)) +...
YI(ii,jj)*sin(VA(ii)-VA(jj)));
QSUM = QSUM + VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj)) -...
YI(ii,jj)*cos(VA(ii)-VA(jj)));
end
PCAL(ii) = PSUM;
QCAL(ii) = QSUM;
end
%End of functionCalculatePowers

-----------------------

function [PCAL,QCAL,ltcPCAL,ltcQCAL] = LTCCalculatedPowers(NLTC,...
LTCsend,LTCrec,Tap,Rltc,Xltc,VM,VA,PCAL,QCAL)
for ii = 1: NLTC
kk = (ii-1)*2+1;
% Calculate LTC admittances
denom = Rltc(ii)^2+Xltc(ii)^2;
YRS = Rltc(ii)/denom;
YIS = -Xltc(ii)/denom;
YRM = -Rltc(ii)/denom;
YIM = Xltc(ii)/denom;
A1 = VA(LTCsend(ii))-VA(LTCrec(ii));
A2 = VA(LTCrec(ii))-VA(LTCsend(ii));
% Calculate LTC powers
ltcPCAL(kk) = VM(LTCsend(ii))^2*YRS + Tap(ii)*VM(LTCsend(ii))*...
VM(LTCrec(ii))*(YRM*cos(A1) + YIM*sin(A1));
ltcQCAL(kk) = -VM(LTCsend(ii))^2*YIS + Tap(ii)*VM(LTCsend(ii))*...
VM(LTCrec(ii))*(YRM*sin(A1) - YIM*cos(A1));
ltcPCAL(kk+1) = (VM(LTCrec(ii))*Tap(ii))^2*YRS + Tap(ii)*...
VM(LTCsend(ii))*VM(LTCrec(ii))*(YRM*cos(A2)+YIM*sin(A2));
ltcQCAL(kk+1) = - (VM(LTCrec(ii))*Tap(ii))^2*YIS + Tap(ii)*...
VM(LTCsend(ii))*VM(LTCrec(ii))*(YRM*sin(A2)-YIM*cos(A2));
% Update calculated powers PCAL and QCAL
PCAL(LTCsend(ii)) = PCAL(LTCsend(ii)) + ltcPCAL(kk);
QCAL(LTCsend(ii)) = QCAL(LTCsend(ii)) + ltcQCAL(kk);
PCAL(LTCrec(ii)) = PCAL(LTCrec(ii)) + ltcPCAL(kk+1);
QCAL(LTCrec(ii)) = QCAL(LTCrec(ii)) + ltcQCAL(kk+1);
end

-------------------------

%Function to compute power mismatches
function [DPQ,DP,DQ,flag] = PowerMismatches(nmax,nbb,tol,bustype,...
flag,PNET,QNET,PCAL,QCAL)
% POWER MISMATCHES
DPQ = zeros(1,nmax);
DP = zeros(1,nbb);
DQ = zeros(1,nbb);
DP = PNET - PCAL;
DQ = QNET - QCAL;
% To remove the active and reactive powers contributions of the slack
% bus and reactive power of all PV buses
for ii = 1: nbb
if (bustype(ii) == 1 )
DP(ii) = 0;
DQ(ii) = 0;
elseif (bustype(ii) == 2 )
DQ(ii) = 0;
end
end
% Re-arrange mismatch entries
kk = 1;
for ii = 1: nbb
DPQ(kk) = DP(ii);
DPQ(kk+1) = DQ(ii);
kk = kk + 2;
end
% Check for convergence
for ii = 1: nbb*2
if ( abs(DPQ) < tol)
flag = 1;
end
end
%End function PowerMismatches

----------------------------------

    %Function to built the Jacobian matrix
    function [JAC] = NewtonRaphsonJacobian(nmax,nbb,bustype,PCAL,QCAL,...
    VM,VA,YR,YI)
    % JACOBIAN FORMATION
    % Include all entries
    JAC = zeros(nmax,nmax);
    iii = 1;
    for ii = 1: nbb
    jjj = 1;
    for jj = 1: nbb
    if ii == jj
    JAC(iii,jjj) = -QCAL(ii) - VM(ii)^2*YI(ii,ii);
    JAC(iii,jjj+1) = PCAL(ii) + VM(ii)^2*YR(ii,ii);
    JAC(iii+1,jjj) = PCAL(ii) - VM(ii)^2*YR(ii,ii);
    JAC(iii+1,jjj+1) = QCAL(ii) - VM(ii)^2*YI(ii,ii);
    else
    JAC(iii,jjj) = VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj))...
    -YI(ii,jj)*cos(VA(ii)-VA(jj)));
    JAC(iii+1,jjj) = -VM(ii)*VM(jj)*(YI(ii,jj)*sin(VA(ii)...
    -VA(jj))+YR(ii,jj)*cos(VA(ii)-VA(jj)));
    JAC(iii,jjj+1) = -JAC(iii+1,jjj);
    JAC(iii+1,jjj+1) = JAC(iii,jjj);
    end
    jjj = jjj + 2;
    end
    iii = iii + 2;
    end
    % Delete the voltage magnitude and phase angle equations of the slack
    % bus and voltage magnitude equations corresponding to PV buses
    for kk = 1: nbb
    if (bustype(kk) == 1)
    ii = kk*2-1;
    for jj = 1: 2*nbb
    if ii == jj
    JAC(ii,ii) = 1;
    else
    JAC(ii,jj) = 0;
    JAC(jj,ii) = 0;
    end
    end
    end
    if (bustype(kk) == 1) || (bustype(kk) == 2)
    ii = kk*2;
    for jj = 1: 2*nbb
    if ii == jj
    JAC(ii,ii) = 1;
    else
    JAC(ii,jj) = 0;
    JAC(jj,ii) = 0;
    end
    end
    end
    end
    %End of function NewtonRaphsonJacobian

---------------------------------

function [JAC] = LTCNewtonRaphsonJacobian(bustype,LTCsend,LTCrec,...
NLTC,Tap,Bus,Rltc,Xltc,ltcPCAL,ltcQCAL,VM,VA,JAC)
% LTC JACOBIAN MODIFICATION
for ii = 1: NLTC
ind = Bus(ii)-LTCsend(ii);
JAC(:,2*Bus(ii))= 0.0;
for nn = 1: 2
% Calculate LTC admittances
denom = Rltc(ii)^2+Xltc(ii)^2;
YRS = Rltc(ii)/denom;
YIS = - Xltc(ii)/denom;
% Calculate LTC Jacobian entries
JKK(1,1) = - (VM(LTCsend(ii))^2)*YIS;
JKK(1,2) = (VM(LTCsend(ii))^2)*YRS;
JKK(2,1) = - (VM(LTCsend(ii))^2)*YRS;
JKK(2,2) = - (VM(LTCsend(ii))^2)*YIS;
JKM(1,1) = ltcQCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YIS;
JKM(1,2) = ltcPCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YRS;
JKM(2,1) = - (ltcPCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YRS);
if ind == 0
JKM(2,2) = -(-ltcQCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YIS);
else
JKM(2,2) = ltcQCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YIS;
end
if ((bustype(LTCsend(ii)) == 4) && (Bus(ii) == LTCsend(ii)) )
JKK(1,2) = (ltcPCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YRS);
if (nn == 2)
JKK(2,2) = -(-ltcQCAL((ii-1)*2+nn)+ (VM(LTCsend(ii))^2)*YIS);
else
JKK(2,2) = ltcQCAL((ii-1)*2+nn) + (VM(LTCsend(ii))^2)*YIS;
JKM(2,1) = - (ltcPCAL((ii-1)*2+nn) + ...
(VM(LTCsend(ii))^2)*YRS);
JKM(2,2) = (ltcQCAL((ii-1)*2+nn) + ...
(VM(LTCsend(ii))^2)*YIS);
end
end
if ( (bustype(LTCsend(ii))==2) && (bustype(LTCrec(ii))>2))
JKK(1,2) = 0;
JKK(2,1) = 0;
JKK(2,2) = 0;
if nn == 1
JKM(2,1) = 0;
JKM(2,2) = 0;
else
JKM(1,2) = 0;
JKM(2,2) = 0;
end
elseif ( (bustype(LTCsend(ii))==1) && (bustype(LTCrec(ii))>2))
JKK = zeros;
JKM = zeros;
JMK = zeros;
end
kk = 2*LTCsend(ii)-1;
mm = 2*LTCrec(ii)-1;
JAC(kk:kk+1,kk:kk+1) = JAC(kk:kk+1,kk:kk+1) + JKK;
JAC(kk:kk+1,mm:mm+1) = JAC(kk:kk+1,mm:mm+1) + JKM;
send = LTCsend(ii);
LTCsend(ii) = LTCrec(ii);
LTCrec(ii) = send;
VM(LTCsend(ii)) = VM(LTCsend(ii))*Tap(ii);
end
end

-----------------------------

%Function to update state variables
function [VA,VM] = StateVariablesUpdates(nbb,D,VA,VM)
iii = 1;
for ii = 1: nbb
VA(ii) = VA(ii) + D(iii);
VM(ii) = VM(ii) + D(iii+1)*VM(ii);
iii = iii + 2;
end
%End function StateVariableUpdating

-----------------------

function [VM,Tap] = LTCUpdates(VM,D,bustype,NLTC,LTCsend,LTCrec,Tap,...
Bus,LTCVM)
for ii = 1: NLTC
if ((bustype(LTCsend(ii)) == 4) && (Bus(ii) == LTCsend(ii)))
Tap(ii) = Tap(ii) + (D(2*LTCsend(ii))*Tap(ii));
VM(LTCsend(ii)) = LTCVM(ii);
elseif ( (bustype(LTCrec(ii)) == 4) && (Bus(ii) == LTCrec(ii)) )
Tap(ii) = Tap(ii) + D(2*LTCrec(ii))*Tap(ii);
VM(LTCrec(ii)) = LTCVM(ii);
end
end

-------------------------

function [Tap,bustype] = LTCLimits(bustype,NLTC,Tap,TapHi,TapLo,...
LTCsend,LTCrec)
% CHECK FOR POSSIBLE LTCs’ TAPS LIMITS VIOLATIONS
for ii = 1: NLTC
for kk = 1: 2
if (bustype(LTCsend(ii)) == 4)
if ( Tap(ii) > TapHi(ii) )
Tap(ii) = TapHi(ii);
bustype(LTCsend(ii)) = 3;
elseif (Tap(ii) < TapLo(ii))
Tap(ii) = TapLo(ii);
bustype(LTCsend(ii)) = 3;
end
end
LTCsend(ii) = LTCrec(ii);
end
end

-------------------------------

function [LTCPQsend,LTCPQrec] = LTCPQflows(NLTC,LTCsend,LTCrec,Rltc,...
Xltc,Tap,VM,VA)
for ii = 1: NLTC
% Calculate LTC admittances
denom = Rltc(ii)^2+Xltc(ii)^2;
YRS = Rltc(ii)/denom;
YIS = -Xltc(ii)/denom;
YRM = -Rltc(ii)/denom;
YIM = Xltc(ii)/denom;
for jj = 1 : 2
A1 = VA(LTCsend(ii))-VA(LTCrec(ii));
% Calculate LTC powers
ltcPCAL = VM(LTCsend(ii))^2*YRS + Tap(ii)*VM(LTCsend(ii))*VM(LTCrec(ii))*(YRM*cos(A1) + YIM*sin(A1));
ltcQCAL = -VM(LTCsend(ii))^2*YIS + Tap(ii)*VM(LTCsend(ii))*VM(LTCrec(ii))*(YRM*sin(A1) - YIM*cos(A1));
if jj == 1
LTCPQsend = ltcPCAL + 1i*ltcQCAL;
else
LTCPQrec = ltcPCAL + 1i*ltcQCAL;
end
send = LTCsend(ii);
LTCsend(ii) = LTCrec(ii);
LTCrec(ii) = send;
end
end

------------------------

%Function to calculate the power flows
function [PQsend,PQrec,PQloss,PQbus] = PQflows(nbb,ngn,ntl,nld,...
genbus,loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,...
QLOAD,VM,VA)
PQsend = zeros(1,ntl);
PQrec = zeros(1,ntl);
% Calculate active and reactive powers at the sending and receiving
% ends of tranmsission lines
for ii = 1: ntl
Vsend = ( VM(tlsend(ii))*cos(VA(tlsend(ii))) + ...
VM(tlsend(ii))*sin(VA(tlsend(ii)))*1i );
Vrec = ( VM(tlrec(ii))*cos(VA(tlrec(ii))) + ...
VM(tlrec(ii))*sin(VA(tlrec(ii)))*1i );
tlimped = tlresis(ii) + tlreac(ii)*1i;
current =(Vsend - Vrec) / tlimped + Vsend*( tlcond(ii) + ...
tlsuscep(ii)*1i )*0.5 ;
PQsend(ii) = Vsend*conj(current);
current =(Vrec - Vsend) / tlimped + Vrec*( tlcond(ii) + ...
tlsuscep(ii)*1i )*0.5 ;
PQrec(ii) = Vrec*conj(current);
PQloss(ii) = PQsend(ii) + PQrec(ii);
end
% Calculate active and reactive powers injections at buses
PQbus = zeros(1,nbb);
for ii = 1: ntl
PQbus(tlsend(ii)) = PQbus(tlsend(ii)) + PQsend(ii);
PQbus(tlrec(ii)) = PQbus(tlrec(ii)) + PQrec(ii);
end
% Make corrections at generator buses, where there is load, in order to
% get correct generators contributions
for ii = 1: nld
jj = loadbus(ii);
for kk = 1: ngn
ll = genbus(kk);
if jj == ll
PQbus(jj) = PQbus(jj) + ( PLOAD(ii) + QLOAD(ii)*1i );
end
end
end
%End function PQflows

-------------------------------

this code is in 'Facts:Modelling and Simulating in Power Networks, by Enrique Acha'. this have warning. 'Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.334107e-038.' I can't fix it. Can anyone help me, please?

[Paulo Silva]

Nice code, too bad it isn't properly formatted so we can test it without having to format it by hand.

[Elif]

i can send you M-files. But i don't know how to send.

[Paulo Silva]

it's good now, no worries, please follow Andrew suggestions

Answer 1

[rida Alturshani]

the problem that you just copy and paste the code this is the right code D=inv(jac)*DPQ'.

Answer 2

[Andrew Newell]

This is a lot of code! You could start by finding where the warning is generated. Type

dbstop if warning

before running the code, or choose the menu option Debug -> Stop if errors/warnings. This will stop you in the program where the warning occurs and you can find out what variables are involved. See debug for how to use the debugger.

[Elif]

i know where the warning occurs but i can't fix it.

Warning from ==> LTCNewtonRaphson at 33

D = JAC\DPQ';

But I didn't know the option you wrote. Thanks for that.

[Andrew Newell]

So you are dividing DPQ by the ill-conditioned matrix JAC. It appears that this matrix is generated by NewtonRaphsonJacobian, which in turn depends on a lot of other things, many of which are themselves calculated by functions. There probably isn't a single source for the warning.

[Andrew Newell]

I think the real problem with this code is that there are an awful lot of parameters. If you set up a linear algebra problem with lots of arbitrary parameters, you're bound to get an ill-conditioned problem some of the time. The warning might go away if you vary some of the parameters.

Answer 3

[Matt Tearle]

The possibility of an ill-conditioned Jacobian is an inherent property of Newton/gradient methods. It occurs whenever the function is locally "flat" (in however many dimensions) at the point you're currently evaluating. There's been a lot of theoretical work to show that the iterates taken by Newton methods are essentially chaotic -- you practically don't know if you're going to hit one of this problem areas.

Courtesy of Paulo Silva, here's some background info on these methods.

[Elif]

I also use newton/raphson method to analyse a system with 5 busses, without tap changing transformer. The code work. I added tap changer as a new bus and new functions. I checked all coded and it seemed correct to me.

Answer 4

[qendrim zeka]

bustype data for these buses should be like this and the code will works, you will get the same results as the book

bustype(3) = 4 ; VM(3) = 1 ; VA(3) =0 ;
bustype(6) = 3; VM(6) = 1 ; VA(6) =0 ;

Answer 5

[Haider Ahmed]

dear elif what the solve for this error?? because I have the same problem.