SetDirectory[NotebookDirectory[]]
SetOptions[EvaluationNotebook[],Background->White]

ClearAll["Global`*"]

funkce

Ymatice[uzel1_,uzel2_,Z_,Y_,B_,T_]:=Module[{Ypom,u1,u2,m},
Ypom=(1/Z);
Ykm=(Ypom/T);
Yk0=(1-(1/T))*Ypom+I*(B/2);
Ym0=((1-T)/(T^2))*Ypom+I*(B/2);
(*Nejprve je nutno vložit na diagnálu uzlové admitance*)
m=DiagonalMatrix[Y];
For[i=1,i<Length[uzel1]+1,i++,
u1=uzel1[[i]];
u2=uzel2[[i]];
(*Dále pak admitance všech okolních větví*)
m[[u1,u1]]=m[[u1,u1]]+Ykm[[i]]+Yk0[[i]];
m[[u2,u2]]=m[[u2,u2]]+Ykm[[i]]+Ym0[[i]];
(*nediagonalni prvky*)
m[[u1,u2]]=m[[u1,u2]]-Ykm[[i]];
m[[u2,u1]]=m[[u2,u1]]-Ykm[[i]];
];
m
];

dos[x_,hod_]:=Table[x[[ii]]->hod[[ii]],{ii,1,Length@x}]

Vstupní data optimalizace

UZEL=Import["uzel.xlsx"][[1]]

UZEL//MatrixForm

(*vztažné napětí*)
Uvzad=UZEL[[2;;Length@UZEL,2]]
(*ve slacku a PU uzlech je to zadáná hodnota, pro PQ je to odhad*)
Uzad=UZEL[[2;;Length@UZEL,3]]
(*hodnoty úhlů*)
δzad=UZEL[[2;;Length@UZEL,4]]
(*PQ-1,PU-2,Slack-3*)
typyUzlu=UZEL[[2;;Length@UZEL,5]]//Round

(*vztažné hodnoty*)
Sv=1;
Zv=((Uvzad^2)/Sv)
Iv=(Sv/((Sqrt[3])*Uvzad))

(*pozice jednotlivých typů*)
PQPoz=Position[typyUzlu,1]/.{x_}->x
PUPoz=Position[typyUzlu,2]/.{x_}->x
SlackPoz=Position[typyUzlu,3]/.{x_}->x

Sv=1;
(*Pdod-Podb*)
PWzad=((UZEL[[2;;Length@UZEL,6]]-UZEL[[2;;Length@UZEL,8]])/Sv);
(*Qdod-Qodb*)
QWzad=((UZEL[[2;;Length@UZEL,7]]-UZEL[[2;;Length@UZEL,9]])/Sv);
(*dosazovací vektor žádaných hodnot výkonů v uzlech PQ a PU*)
Table[(Subscript[P, ii])->%%[[ii]],{ii,Join[PQPoz,PUPoz]}];
PQWDos=Join[%,Table[(Subscript[Q, ii])->%%[[ii]],{ii,PQPoz}]]

(*Pevné hodnoty napětí úhlu δ, které nechceme zjišťovat.*)
(*pro slack je to úhel δ*)
δUPevnaDos={Table[{(Subscript[δ, uzel])->δzad[[uzel]]},{uzel,SlackPoz}]
}//Flatten

(*promenne*)
Join[
(*pro všechny uzly typu PQ a PU hledáme úhel*)
Table[(Subscript[δ, uzel]),{uzel,Join[PQPoz,PUPoz]}],
(*velikosti napětí hledáme pouze pro PQ uzly*)
Table[(Subscript[U, uzel]),{uzel,PQPoz}]
];
δU=%
PQ=%/.{δ->P,U->Q}

UU=Table[(Subscript[U, i])*Exp[I*(Subscript[δ, i])],{i,1,Length@Uzad}]

(*počáteční odhady*)
ξ0Dos={
Table[(Subscript[δ, uzel])->δzad[[uzel]],{uzel,Join[PQPoz,PUPoz]//Sort}],
Table[(Subscript[U, uzel])->Uzad[[uzel]],{uzel,1,Length@Uzad}],
Table[(Subscript[λP, ii])->0,{ii,Join[{R},PQPoz]}],
Table[(Subscript[λQ, ii])->0.1,{ii,PQPoz}]
}//Flatten
%//Length

ξ={
Table[(Subscript[δ, uzel]),{uzel,Join[PQPoz,PUPoz]//Sort}],
Table[(Subscript[U, uzel]),{uzel,1,Length@Uzad}],
Table[(Subscript[λP, ii]),{ii,Join[{R},PQPoz]}],
Table[(Subscript[λQ, ii]),{ii,PQPoz}]
}//Flatten
%//Length

Popis sítě - admitanční matice

VETEV=Import["vetev.xlsx"][[1]];
%//MatrixForm

(*propoje mezi uzly*)
uzel1=VETEV[[2;;Length@VETEV,1]]//Round
uzel2=VETEV[[2;;Length@VETEV,2]]//Round

(*vektor uzlových admitancí sousledné složky*)
Y1=ConstantArray[0,Max[Join[uzel1,uzel2]]]

(*vektor meziuzlových impedancí a susceptancí sousledné složky*)
R1=VETEV[[2;;Length@VETEV,3]];
X1=VETEV[[2;;Length@VETEV,4]];
Zkm1=((R1+I*X1)/(Zv[[uzel1]]))
B1=VETEV[[2;;Length@VETEV,5]]*Zv[[uzel1]]*(10^(-6))
T=VETEV[[2;;Length@VETEV,6]]

ymat=Ymatice[uzel1,
uzel2,
Zkm1,
Y1,
B1,
T];
%//MatrixForm

Vyjádření výkonů

SS=DiagonalMatrix[UU].Conjugate[ymat].Conjugate[UU];
%//MatrixForm

SSDos=Table[{(Subscript[P, uzel])->ComplexExpand[Re[SS[[uzel]]]],(Subscript[Q, uzel])->ComplexExpand[Im[SS[[uzel]]]]},{uzel,1,Length@UU}]//Flatten;

Tvorba Lagrangerovy funkce

(Subscript[λP, R])*(Table[(((Subscript[P, ii])(*/.SSDos*))-((Subscript[P, ii])/.PQWDos)),{ii,Join[PUPoz]//Sort}]/.List->Plus)
%+Table[(Subscript[λP, ii])*(((Subscript[P, ii])(*/.SSDos*))-((Subscript[P, ii])/.PQWDos)),{ii,Join[PQPoz]//Sort}]/.List->Plus
PP=%/.SSDos;
%//Length

Table[(Subscript[λQ, ii])*(((Subscript[Q, ii])(*/.SSDos*))-((Subscript[Q, ii])/.PQWDos)),{ii,PQPoz}]/.List->Plus
QQ=%/.SSDos;
%//Length

WD=2000000;
WH=2000000;
Plot[{WD*((Max[0,-(Subscript[U, ii])+0.95])^2),WH*((Max[0,(Subscript[U, ii])-1.05])^2)},{(Subscript[U, ii]),0,1.5},PlotRange->{{0,1.2},{0,1000}},PlotLegends->{"min","max"}]

funcUU=Table[{((WD*(Subscript[lD, ii]))/2)*(((-(Subscript[U, ii])+0.95))^2),((WH*(Subscript[lH, ii]))/2)*((((Subscript[U, ii])-1.05))^2)},{ii,1,Length@Uzad}]//Flatten//Total

WS=200000;
pridanaStabilizace=Table[(1/2)*(WS/(10^iterace))*((((Subscript[U, ii])-1))^2),{ii,Join[PUPoz,PQPoz,SlackPoz]}]/.List->Plus

První iterace - 0

omezUU=Table[{If[(Subscript[U, ii])≥0.95,(Subscript[lD, ii])->0,(Subscript[lD, ii])->1]/.ξ0Dos,If[(Subscript[U, ii])≤1.05,(Subscript[lH, ii])->0,(Subscript[lH, ii])->1]/.ξ0Dos},{ii,1,Length@Uzad}]//Flatten

(*Lagrangerova funke*)
Lobec=((Subscript[P, 1])/.SSDos)+PP+QQ+funcUU+pridanaStabilizace;

(*Gradient*)
GLobec=D[Lobec,{ξ}];
GL0=GLobec/.ξ0Dos/.δUPevnaDos/.omezUU/.iterace->0
%//Length
TableForm[%%,TableHeadings->{Table[ii,{ii,1,Length@ξ}],None}];

(*Ovroubený Hessián*)
Wobec=D[D[Lobec,{ξ}],{ξ}];
W0=Wobec/.ξ0Dos/.δUPevnaDos/.omezUU/.iterace->0;
IW0=%//Inverse;
TableForm[%,TableHeadings->{Table[ii,{ii,1,58}],Table[ii,{ii,1,59}]}];

(*struktůra ovroubeného hesiánu*)
(*
...
...
*)

W0[[1;;31,1;;31]];(*hesova matice*)
W0[[32;;58,1;;31]]==Transpose[W0[[1;;31,32;;58]]]
W0[[32;;58,32;;58]];

Δξ0=-Inverse[W0].GL0

(ξ/.ξ0Dos)+Δξ0

ξkDos=dos[ξ,(ξ/.ξ0Dos)+Δξ0]

Další iterace - k

omezUU=Table[{If[(Subscript[U, ii])≥0.95,(Subscript[lD, ii])->0,(Subscript[lD, ii])->1]/.ξkDos,If[(Subscript[U, ii])≤1.05,(Subscript[lH, ii])->0,(Subscript[lH, ii])->1]/.ξkDos},{ii,1,Length@Uzad}]//Flatten

(*Gradient*)
GLk=GLobec/.ξkDos/.δUPevnaDos/.omezUU
%//Length
TableForm[%%,TableHeadings->{Table[ii,{ii,1,Length@ξ}],None}];

(*Ovroubený Hessián*)
Wk=Wobec/.ξkDos/.δUPevnaDos/.omezUU/.iterace->1;
IWk=%//Inverse;
TableForm[%,TableHeadings->{Table[ii,{ii,1,59}],Table[ii,{ii,1,59}]}];

Δξk=-Inverse[Wk].GLk

ξkDos=dos[ξ,(ξ/.ξkDos)+Δξk]

newton[{k_,odhad_,defekt_}]:=Block[{omezUU,GLk,Wk,Δξk,ξkDos=odhad,iterace=k},
omezUU=Table[{If[(Subscript[U, ii])≥0.95,(Subscript[lD, ii])->0,(Subscript[lD, ii])->1]/.ξkDos,If[(Subscript[U, ii])≤1.05,(Subscript[lH, ii])->0,(Subscript[lH, ii])->1]/.ξkDos},{ii,1,Length@Uzad}]//Flatten;
GLk=GLobec/.ξkDos/.δUPevnaDos/.omezUU;
Wk=Wobec/.ξkDos/.δUPevnaDos/.omezUU;
Δξk=(*ConstantArray[0,Length@ξ]*)-Inverse[Wk].GLk;
ξkDos=dos[ξ,(ξ/.ξkDos)+Δξk];
{k+1,ξkDos,Max[GLk]}
]

NestWhileList[newton,{0,ξ0Dos,1000},Max[#[[3]]]>(10^(-8))&];
MatrixForm[%,TableHeadings->{None,{"k","ξkDos","Max[GLk]"}}]
vysDos=%%[[Length@%%,2]];

(SS/.vysDos/.δUPevnaDos)[[12]]-PWzad[[12]]
(SS/.vysDos/.δUPevnaDos)[[14]]-PWzad[[14]]

SS/.vysDos/.δUPevnaDos
%[[12]]

SS/.vysDos/.δUPevnaDos/.List->Plus

(*ztráty*)
(*
...
...
*)