# Simple 4 bus Newton-Raphson matlab not converging

Tags:
1. Oct 27, 2015

### tim9000

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

%{
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. Oct 27, 2015

### BvU

Hi,

How about putting $[$Code = matlab $]$ and $[$/Code$]$ around your source text ?

3. Oct 27, 2015

### tim9000

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

4. Oct 28, 2015

### BvU

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

%{
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 ?

5. Oct 28, 2015

### tim9000

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!

6. Oct 28, 2015

### BvU

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

7. Oct 28, 2015

### tim9000

No, didn't work, but you still might have made a good point, cheers.

8. Oct 28, 2015

### tim9000

(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
9. Oct 29, 2015

### BvU

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 ?

10. Oct 29, 2015

### tim9000

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.

11. Oct 29, 2015

### BvU

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 ?

12. Oct 29, 2015

### tim9000

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
13. Oct 29, 2015

### BvU

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 .

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$$
 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 ?

--

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
14. Oct 29, 2015

### tim9000

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))) ?

15. Oct 30, 2015

### BvU

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.

16. Oct 30, 2015

### tim9000

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

17. Oct 30, 2015

### BvU

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}$$
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 ?

18. Oct 30, 2015

### tim9000

But when I ran it I got:

19. Oct 30, 2015

### BvU

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^*$

20. Oct 30, 2015

### BvU

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 ?