1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Simple 4 bus Newton-Raphson matlab not converging

  1. Oct 27, 2015 #1
    1. The problem statement, all variables and given/known data

    All the necessary data is in the code, I'm just trying to converge NR, I decided to use the equation S = V^2 / Z since I had the admittance matrix and powers (needed voltages)

    I think my simple algorithm has a slight issue I can't find.

    2. Relevant equations
    Thank you!

    3. The attempt at a solution


    % ** ENTER THE AMOUNT OF BUS NODES FOR THE SIZE OF THE NETWORK **
    bus_num = 4 % FOR A N = 4 BUS NETWORK
    initial_guess = 1 % Guess voltage in P.U.
    slack_bus = 1 % Which bus number is the slack bus
    convergence_criterion = input('What do you want the convergence criteria to be? (difference < 1): '); %assignment spec was 0.0001
    Imax = input('What do you want the maximum iteration to be? ');
    volts = [1; 1; 1; 1]; % Known P.U. bus voltages (inc slack bus Estimation)
    real_powers = [0; 0.2; 0.02; 0.2]; % Known P.U. bus P
    reactive_powers = [0; 0.15; 0.15; 0.15];% Known P.U. bus Q
    Conj_power = conj(complex(real_powers, reactive_powers)); % *** Add Power to create Conjugate Apparent Power ***

    bus_row1 = [0, 0.5882-j*2.3529, 0, 1.1765-j*4.7059]; % ENTER THE
    bus_row2 = [0.5882-j*2.3529, 0, 0.3846-j*1.9231, 0.5882-j*2.3529]; % ADMITTANCES
    bus_row3 = [0, 0.3846-j*1.9231, 0, 1.1765-j*4.7059]; % BETWEEN BUSES
    bus_row4 = [1.1765-j*4.7059, 0.5882-j*2.3529, 1.1765-j*4.7059, 0]; % FOR MATRIX BELOW

    %{
    The above variables represent admittances
    corresponding to the following matrix:
    [Y11 Y12 Y13 Y14]
    |Y21 Y22 Y23 Y24|
    Y_mat = |Y31 Y32 Y33 Y34|
    [Y41 Y42 Y43 Y44]
    User will need define more 'bus_row' variables
    as above and add them to the Y_mat below if there
    is more than 4 buses:
    %}
    Y_mat = [bus_row1; bus_row2; bus_row3; bus_row4];

    % *** Make the appropreate Admittance matrix for the network: ***
    for r = 1 : bus_num,
    for c = 1 : bus_num,
    if (r ~= c)
    Y_mat(r,r) = Y_mat(r,r) + Y_mat(r,c);
    end
    end
    end
    % *** Turn the Admittances that arn't on
    % the diagonal negative: ***
    for I = 1 : bus_num,
    for k = 1 : bus_num,
    if (I ~= k)
    Y_mat(I,k) = -1*Y_mat(I,k);
    end
    end
    end

    % ** Newton-Raphson Method ***
    volts = [1; 1; 1; 1]; % P.U. bus voltages (inc slack bus Estimation)
    last_guess = [0; 0; 0; 0];
    itterations = 0;
    tol = max(abs(volts - last_guess));


    % (volts)^2)*Y_mat = Conj_power(I)
    % so:
    % function = (volts)^2)*Y_mat
    % constant is Conj_power(I)
    % Jacobian is (2*volts*Y_mat)
    while (tol > convergence_criterion) & itterations < Imax

    last_guess = volts; % FOR TESTING CONVERGENCE

    for I = 1 : bus_num,
    if (I ~= slack_bus)
    for k = 1 : bus_num,
    C3 = 0;
    if (I ~= k)

    C3 = (Conj_power(I) - ((volts(k))^2)*Y_mat(I,k))/(2*volts(k)*Y_mat(I,k));


    end
    end
    volts(I) = last_guess(I) + C3; % Residual I believe
    tol = max(abs(volts - last_guess));
    end
    end

    itterations = itterations + 1;
    end
    % *** Display converged result ***
    if itterations >= Imax,
    disp('Did not converge Newton-Raphson within itteration limit');
    else
    fprintf('\nBonus: The number of itterations required for Newton-Raphson was: %d.\n\n', itterations)
    for I = 1 : bus_num,
    fprintf('The voltage for bus %d is %0.4f V real and %0.4f V reactive, with angle %0.2f degrees \n', I, real(volts2(I, 1)), imag(volts2(I, 1)), angle(volts2(I, 1))*180/pi)
    if (angle(volts2(I, 1))*180/pi >= 45),
    disp('\n This solution seems dodgy, might be wrong convergence answer, try a new initial guess sunshine!')
    end
    end
    end
     
  2. jcsd
  3. Oct 27, 2015 #2

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Hi,

    How about putting ##[##Code = matlab ##]## and ##[##/Code##]## around your source text ?
     
  4. Oct 27, 2015 #3
    Hi,
    Sorry, I should have, I'ven't ever posted code on here. It's too late for me to edit the post now, but I did say it in the title of the thread.
    Cheers
     
  5. Oct 28, 2015 #4

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Code (Matlab M):

    % ** ENTER THE AMOUNT OF BUS NODES FOR THE SIZE OF THE NETWORK **
    bus_num = 4 % FOR A N = 4 BUS NETWORK
    initial_guess = 1 % Guess voltage in P.U.
    slack_bus = 1 % Which bus number is the slack bus
    convergence_criterion = input('What do you want the convergence criteria to be? (difference < 1): '); %assignment spec was 0.0001
    Imax = input('What do you want the maximum iteration to be? ');
    volts = [1; 1; 1; 1]; % Known P.U. bus voltages (inc slack bus Estimation)
    real_powers = [0; 0.2; 0.02; 0.2]; % Known P.U. bus P
    reactive_powers = [0; 0.15; 0.15; 0.15];% Known P.U. bus Q
    Conj_power = conj(complex(real_powers, reactive_powers)); % *** Add Power to create Conjugate Apparent Power ***

    bus_row1 = [0, 0.5882-j*2.3529, 0, 1.1765-j*4.7059]; % ENTER THE
    bus_row2 = [0.5882-j*2.3529, 0, 0.3846-j*1.9231, 0.5882-j*2.3529]; % ADMITTANCES
    bus_row3 = [0, 0.3846-j*1.9231, 0, 1.1765-j*4.7059]; % BETWEEN BUSES
    bus_row4 = [1.1765-j*4.7059, 0.5882-j*2.3529, 1.1765-j*4.7059, 0]; % FOR MATRIX BELOW

    %{
    The above variables represent admittances
    corresponding to the following matrix:
    [Y11 Y12 Y13 Y14]
    |Y21 Y22 Y23 Y24|
    Y_mat = |Y31 Y32 Y33 Y34|
    [Y41 Y42 Y43 Y44]
    User will need define more 'bus_row' variables
    as above and add them to the Y_mat below if there
    is more than 4 buses:
    %}
    Y_mat = [bus_row1; bus_row2; bus_row3; bus_row4];

    % *** Make the appropreate Admittance matrix for the network: ***
    for r = 1 : bus_num,
    for c = 1 : bus_num,
    if (r ~= c)
    Y_mat(r,r) = Y_mat(r,r) + Y_mat(r,c);
    end
    end
    end
    % *** Turn the Admittances that arn't on
    % the diagonal negative: ***
    for I = 1 : bus_num,
    for k = 1 : bus_num,
    if (I ~= k)
    Y_mat(I,k) = -1*Y_mat(I,k);
    end
    end
    end

    % ** Newton-Raphson Method ***
    volts = [1; 1; 1; 1]; % P.U. bus voltages (inc slack bus Estimation)
    last_guess = [0; 0; 0; 0];
    itterations = 0;
    tol = max(abs(volts - last_guess));


    % (volts)^2)*Y_mat = Conj_power(I)
    % so:
    % function = (volts)^2)*Y_mat
    % constant is Conj_power(I)
    % Jacobian is (2*volts*Y_mat)
    while (tol > convergence_criterion) & itterations < Imax

    last_guess = volts; % FOR TESTING CONVERGENCE

    for I = 1 : bus_num,
    if (I ~= slack_bus)
    for k = 1 : bus_num,
    C3 = 0;
    if (I ~= k)

    C3 = (Conj_power(I) - ((volts(k))^2)*Y_mat(I,k))/(2*volts(k)*Y_mat(I,k));


    end
    end
    volts(I) = last_guess(I) + C3; % Residual I believe
    tol = max(abs(volts - last_guess));
    end
    end

    itterations = itterations + 1;
    end
    % *** Display converged result ***
    if itterations >= Imax,
    disp('Did not converge Newton-Raphson within itteration limit');
    else
    fprintf('\nBonus: The number of itterations required for Newton-Raphson was: %d.\n\n', itterations)
    for I = 1 : bus_num,
    fprintf('The voltage for bus %d is %0.4f V real and %0.4f V reactive, with angle %0.2f degrees \n', I, real(volts2(I, 1)), imag(volts2(I, 1)), angle(volts2(I, 1))*180/pi)
    if (angle(volts2(I, 1))*180/pi >= 45),
    disp('\n This solution seems dodgy, might be wrong convergence answer, try a new initial guess sunshine!')
    end
    end
    end
     
    Makes a difference eh ?
     
  6. Oct 28, 2015 #5
    Indeed!
    I'd no idea it could do that, pitty it doesn't keep the tabs though.

    So did you have any opinion on why it's not converging within the limit?
    Or any issue with my:
    % (volts)^2)*Y_mat = Conj_power(I)
    % so:
    % function = (volts)^2)*Y_mat
    % constant is Conj_power(I)
    % Jacobian is (2*volts*Y_mat)


    C3 = (Conj_power(I) - ((volts(k))^2)*Y_mat(I,k))/(2*volts(k)*Y_mat(I,k));


    ?
    Cheers!
     
  7. Oct 28, 2015 #6

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    What about the location of C3 = 0 ? Do you want it inside the loop ?

    Note: I'm looking at this as an "interested party". No idea what this is about - but I can guess a little. And, alas, also no MatLab expert
     
  8. Oct 28, 2015 #7
    No, didn't work, but you still might have made a good point, cheers.
     
  9. Oct 28, 2015 #8
    (after moving C3 outside the k loop) I think one thing that might still be wrong is my jacobian.
    From what I remember NR uses the form y = f(x)
    where y is a constant (in this case the powers) and x is the voltages ('volts' matrix), so my function is
    (volts)^2)*Y_mat

    so the jabobian is:
    (2*volts*Y_mat)

    ? Or have I not differentiated that properly?

    I'd like to use Complex power = complex voltage * complex admittance
    like P + Q = V^2 / (R + X)
    but I'm not sure if that relation is true just for apparent power and absolute voltages and impedances?
     
    Last edited: Oct 28, 2015
  10. Oct 29, 2015 #9

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    NR does ##\vec x_{n+1} = \vec x_n + ( \vec {\vec J_n})^{-1} \, \vec r_n## where ##\vec r_n## is the residual vector. Can you distinguish those in your routine ?
     
  11. Oct 29, 2015 #10
    Hey BvU,
    yeah sure:
    so volts is \vec x_{n+1} (sorry the maths scrip I copied from you didn't work, volts is xn+1)
    last_guess is \vec x_n (xn)
    the residual vector is something like:
    conj(Conj_power(I)) - ((volts(k))^2)*Y_mat(I,k)

    and the Jacobian is: (2*volts*Y_mat)

    But I'm not confindent in those, especially the jacobian, and if the power has to be absolute or not. In what I just said then, I un-conjugated the power to be the original power.
    When I did gauss-seidal I had the power negative because they were loads.
     
  12. Oct 29, 2015 #11

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    You use Latex by writing ##\#\# ## in front and closing off with ##\#\# ## :

    ##\#\# ## \vec x_{n+1} ##\#\# ## thus becomes ## \vec x_{n+1} ##

    Alternatively you can open and close with ##$$## to get centered displayed equations (more room for fractions and such);
    ##$$## \vec x_{n+1} ##$$ ## thus becomes $$ \vec x_{n+1} $$

    With ##(\vec {\vec J}_n)^{-1}## I mean the inverse of the Jacobian. I don't see you inverting anywhere ?

    For NR, you normally first evaluate the residual and then apply the inverse of the Jacobian. You can't do it in one loop !

    Also: Volts is a vector, but C3 ?
     
  13. Oct 29, 2015 #12
    Thanks for the Latex pointer.

    I take your point, but doesn't C3 only need to be one complex value because it's going into volts(I) ? (one element)
    Similarly, isn't the derivative: 2*volts*Y_mat
    only one element big, so it can't be inverted?

    Hmm, Is there a better way to do this, rather than evaluating one element at a time?
    It's just that I'm confused because the voltage and powers are both 4x1 matrixes, whereas the admittances are 4x4.

    Anyway, isn't dividing by it the same as multiplying by the inverse?
     
    Last edited: Oct 29, 2015
  14. Oct 29, 2015 #13

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    To help you further, I suppose it would be relevant to know what this is about. I found this (JD McCalley, Iowa State) with a lot of stuff like yours, but for your case I still can't spell out what exactly is being solved for. In short: the problem statement, all variables and given/known data and the relevant equations from the template :smile:.

    I see you putting together matrix Y somewhat as in eq (6) in the link, but you keep the diagonal elements at zero. Why ?

    His (9) ##\vec I = \vec{\vec Y} \vec V## and the power vector is (10), ##\ \ S_k = V_k I_k^*\ ##;

    ## \vec S^*## could well be your conj_power which is a known vector of 4 complex values.
    (why take the conjugate ? It complicates unnecessarily ?).​

    I get the impression you want to solve his eq (12) using NR: $$ \vec S - \vec V \left ( {\vec{\vec Y}} \,\vec V \right )^* = \vec 0 $$
    [edit] I had an inner product ##\vec V\cdot \vec I^*## here, but it should have been ##\ \ S_k = V_k I_k^*\ ## -- so I removed the ## \cdot ## (the dot).​


    NR is about solving N equations in N unknowns, with N = 4 here,
    written in the form ##\vec{\vec f} (\vec x) = \vec 0##

    NR starts with an intial guess ##{\vec x}_0 ## and calculates an initial residual ##{\vec r}_0##.

    and iterates according to ##{\vec x}_{n+1} = {\vec x}_{n} - ({\vec{\vec J}}_n)^{-1} \, \vec r_n ##

    where ## {\vec{\vec J} }## is the Jacobian -- ## J_{ij} = {\partial f_i \over \partial x_j} ## and J sure isn't a number !

    Also for you, all variables and f elements are complex, so you need the complex derivative ! And still run the risk of the associated convergence problems.

    Perhaps you want to consider using his (15) instead ?

    --

    [edit]
    I missed that C3 is used four times, yes. But first you need a residue vector, then you can apply the inverse of the Jacobian on it. At least, that's what I'm used to.
     
    Last edited: Oct 30, 2015
  15. Oct 29, 2015 #14
    Thanks for your continued support.

    I had it as a conjugate for when I had to work out the Gauss-seidal, yeah it is an extra complication, but I just keep it in mind.

    Are you sure? In't it:
    My final Y_mat =
    1.7647 - 7.0588i -0.5882 + 2.3529i 0.0000 + 0.0000i -1.1765 + 4.7059i
    -0.5882 + 2.3529i 1.5610 - 6.6289i -0.3846 + 1.9231i -0.5882 + 2.3529i
    0.0000 + 0.0000i -0.3846 + 1.9231i 1.5611 - 6.6290i -1.1765 + 4.7059i
    -1.1765 + 4.7059i -0.5882 + 2.3529i -1.1765 + 4.7059i 2.9412 -11.7647i

    Because I used these: [Code = matlab]
    [/Code]
    % *** Make the appropreate Admittance matrix for the network: ***
    for r = 1 : bus_num,
    for c = 1 : bus_num,
    if (r ~= c)
    Y_mat(r,r) = Y_mat(r,r) + Y_mat(r,c);
    end
    end
    end
    % *** Turn the Admittances that arn't on
    % the diagonal negative: ***
    for I = 1 : bus_num,
    for k = 1 : bus_num,
    if (I ~= k)
    Y_mat(I,k) = -1*Y_mat(I,k);
    end
    end
    end
    [Code = matlab]
    [/Code]

    I'm trying to find the bus voltages (the right complex values for the volts matrix) hence why I treated it as the variable.
    real_powers = [0; 0.2; 0.02; 0.2]; % Known P.U. bus P
    reactive_powers = [0; 0.15; 0.15; 0.15];% Known P.U. bus Q
    You have seen the admittance matrix.

    Indeed.
    I've tried:
    As:
    C3 = (conj(Conj_power(I)) - ((volts(k)))*conj(volts(k)*(Y_mat(I,k)))) / conj((2*(volts(k))*(Y_mat(I,k))));

    and I've also tried S - V2 / Z*
    as:
    C3 = (conj(Conj_power(I)) - ((volts(k))^2)*conj(Y_mat(I,k))) / ((2*(volts(k))*conj(Y_mat(I,k))));

    I've also tried making the denominator (prime of function) negative.
    But surely if the function is volts(I)^2 then the derivitive of that is 2*volts(I) ?
    But I alreay know the powers?

    But isn't my residual vector just like (conj(Conj_power(I)) - ((volts(k))^2)*conj(Y_mat(I,k))) ?

    Thanks for your assistance!!
     
  16. Oct 30, 2015 #15

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    But surely if the function is volts(I)^2 then the derivitive of that is 2*volts(I) ?

    But the function is not so simple. There is a whole matrix Y in there. See his eq (12)​

    ##[##\code##]## should be at the end of the code section and = matlab probably without spaces

    But I alreay know the powers?

    I mean: use his (15) as the equations to solve for

    Sory, short on time; have to run.

    :smile:
     
  17. Oct 30, 2015 #16
    Hey.

    Oh boy this is driving me crazy.

    Well I can interpret his (12) to make my function:

    (conj(Conj_power(I)) - (volts(I)*conj(volts(k))*conj(Y_mat(I,k)))) = 0

    where I is the outer loop and K is the inner loop. But I don't know what the Jacobian would be for that?

    I'm unfortunately still at a loss as to how I could use his equation (15) to solve for the voltages though.

    Thanks
     
  18. Oct 30, 2015 #17

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    After Y_mat = [bus_row1; bus_row2; bus_row3; bus_row4], doesn't Y look like $$
    \begin{pmatrix}
    0 & 0.5882-j*2.3529 & 0 & 1.1765-j*4.7059\\
    0.5882-j*2.3529 & 0 & 0.3846-j*1.9231 & 0.5882-j*2.3529\\
    0 & 0.3846-j*1.9231 & 0 & 1.1765-j*4.7059\\
    1.1765-j*4.7059 & 0.5882-j*2.3529 & 1.1765-j*4.7059 & 0
    \end{pmatrix} $$
    [edit]don't understand. Latex preview looks reasonable, browser doesn't do the TeX for some reason.
    Anyway, diagonal elements are zero. Then you have two loops, one that skips r=c and one that skips I=k, so how can the diagonal elements change ?
     
  19. Oct 30, 2015 #18

    But when I ran it I got:
     
  20. Oct 30, 2015 #19

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Then: the residue calculation is still not right. There four Conj_power (your S conjugated ?). You want to find four elements of C3
    $$
    \vec C3_i = \vec S_i - \vec V_i \,\left ( {\vec{\vec Y}} \,\vec V \right )^*_i =
    $$to find ##\left ( {\vec{\vec Y}} \,\vec V \right )^*_i## your inner loop is over j: ##\displaystyle \sum_j Y_{ij}^* V_j^*##
     
  21. Oct 30, 2015 #20

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Yes, you posted that in #14 too. But that doesn't make me understand it. So you have to help me there. How can the diagonal elements be nonzero ?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Simple 4 bus Newton-Raphson matlab not converging
  1. Newton-Raphson help! (Replies: 7)

Loading...