close all
clear all
clc
P= 90;
U1= 220;
U2= 230;
Q1= 30;
Q2= 35;
R13= 3;
R23= 4;
eps= 0.005;
xi=0.1;
b1a=0.005;
b1b=0.25;
b1c=4;
b2a=0.003;
b2b=0.12;
b2c=3;

A= [(b1a-b2a) (b1b+b2b+b2a*2*P) (b1c-b2c-b2b*P-b2a*P^2)];
Aa= roots(A)
P1_0= 34.809 %tipovana hodnota (ph=+ tak zvysujem) (ph=- tak znizujem) %vysledok 35,0750
P2_0= P-P1_0

% I. iteracia
Pb_0= P2_0;
Qb= Q2;
Ub= U2;
deltaP_0= ((P1_0^2+Q1^2)/U1^2)*R13 + ((Pb_0^2+Qb^2)/Ub^2)*R23

W_0= P1_0 + Pb_0 - P - deltaP_0 % bilancia činných výkonov prva presnost epsilon nesplnena

deltaPb_0= -(W_0/(1-(2*Pb_0*R23/(U2^2))))
Pb_1= Pb_0 + deltaPb_0
deltaP_1= ((P1_0^2+Q1^2)/U1^2)*R13 + ((Pb_1^2+Qb^2)/Ub^2)*R23
W_1= P1_0 + Pb_1 - P - deltaP_1 % druha presnost epsilon splnena

deltaPb_1= -(W_1/(1-(2*Pb_1*R23/(U2^2))))
Pb_2= Pb_1 + deltaPb_1
deltaP_2= ((P1_0^2+Q1^2)/U1^2)*R13 + ((Pb_2^2+Qb^2)/Ub^2)*R23
W_2= P1_0 + Pb_2 - P - deltaP_2 % druha presnost epsilon splnena

b1= 4 + 0.25*P1_0 + 0.005*P1_0^2 %menim B1C, B1B, B1A
b2= 3 + 0.12*Pb_2 + 0.003*Pb_2^2 %menim B2C, B2B, B2A 

p1= ((2*P1_0*R13/(U1^2))-(2*Pb_2*R23/(Ub^2)))/(1-(2*Pb_2*R23/(Ub^2)))

bb= b2;
ph= bb - (b1/(1-p1)) % podmienka hospodarnosti splnena

b1n= bb*(1-p1)

B= [0.005  0.25  4-b1n]; %meniš podla b1 riadku B1A, B1B, B1C
Bb= roots(B);
P1_1= Bb(2,1)

P2_1= P+deltaP_1-P1_1
%tipujem podla podmienky hospodarnosti vo výsledkoch ph a sledujem ci je
%kladná alebo záporná