Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
mf=MatrixForm;
randomMatrix[rnm_]:=Table[RandomReal[{-1,1}]&[i],{i,rnm},{j,rnm}];

(*http://kfe.fjfi.cvut.cz/~limpouch/numet/nr/node14.html*)
(*http://kfe.fjfi.cvut.cz/~limpouch/numet/nr/node15.html*)

(*load flow problem, cartesian form of complex numbers*)
Quiet@Remove["Global`*"];
simplify=#&;(*"lower" programing languages have no Mathematica Simplify[#]&*)
(*simplify=Simplify;*)
conjugate[{re_,im_}]:={re,-im};
times[ua:{uare_,uaim_},ub:{ubre_,ubim_}]:={uare*ubre-uaim*ubim,uare*ubim+ubre*uaim};
re[{re_,im_}]:={re,0};
im[{re_,im_}]:={im,0};
sqabs[{re_,im_}]:={(re^2)+(im^2),0};
powerMinus1[z:{re_,im_}]:=(({re,-im})/((re^2)+(im^2)));
ratio[a_,b_]:=times[powerMinus1[b],a]
apparentPower[u_,i_]:=times[u,conjugate[i]];
activePower[u_,i_]:=re[apparentPower[u,i]];
reactivePower[u_,i_]:=im[apparentPower[u,i]];
current[u1_,u2_,z_]:=ratio[u1-u2,z];
toMyForm[z_]:={Re[z],Im[z]};
fromMyForm[{re_,im_}]:=re+I*im;
rearrangeRealEqn[eqn_Equal]:={simplify[eqn]}/.Equal[a_,b_]:>Equal[a-b,0];
rearrangeComplexEqn[{rea_,ima_}=={reb_,imb_}]:=DeleteCases[#,rearrangeRealEqn[True]]&@Flatten[rearrangeRealEqn/@{simplify[rea==reb],simplify[ima==imb]}];
rearrangeComplexEqn[{a,b}=={c,d}];
rearrangeRealEqn[x+2y==y-x];

{(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]),(Subscript[z, 4])}=kz*(toMyForm/@{.01+0.05I,2*(01+0.05I),0.01+0.02I,0.01+0.02I});
{(Subscript[u, #])={(Subscript[ure, #]),(Subscript[uim, #])},(Subscript[i, #])={(Subscript[ire, #]),(Subscript[iim, #])}}&/@Range[6];
powerFactor2=0.8;
powerFactor3=0.7;
S1=(10^3)*{powerFactor2,Sin[ArcCos[powerFactor2]]};
S2=500*{powerFactor3,Sin[ArcCos[powerFactor3]]};
Unom=230.;
eqnZ1=rearrangeComplexEqn[(Subscript[i, 1])==current[(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[z, 1])]];
eqnZ2=rearrangeComplexEqn[(Subscript[i, 2])==current[(Subscript[u, 2]),(Subscript[u, 3]),(Subscript[z, 2])]];
eqnZ3=rearrangeComplexEqn[(Subscript[i, 3])==current[(Subscript[u, 5]),(Subscript[u, 3]),(Subscript[z, 3])]];
eqnZ4=rearrangeComplexEqn[(Subscript[i, 4])==current[(Subscript[u, 4]),(Subscript[u, 3]),(Subscript[z, 4])]];
eqnNod1=rearrangeComplexEqn[(Subscript[i, 1])==(Subscript[i, 5])+(Subscript[i, 2])];
eqnNod2=rearrangeComplexEqn[(Subscript[i, 2])+(Subscript[i, 3])+(Subscript[i, 4])==(Subscript[i, 6])];
eqnApparentPower1=rearrangeComplexEqn[apparentPower[(Subscript[u, 2]),(Subscript[i, 5])]==S1];
eqnApparentPower2=rearrangeComplexEqn[apparentPower[(Subscript[u, 3]),(Subscript[i, 6])]==S2];
eqnSourceActivePowerRatio1=rearrangeComplexEqn[activePower[(Subscript[u, 1]),(Subscript[i, 1])]==2*activePower[(Subscript[u, 4]),(Subscript[i, 4])]];
eqnSourceActivePowerRatio2=rearrangeComplexEqn[activePower[(Subscript[u, 1]),(Subscript[i, 1])]==1.5*activePower[(Subscript[u, 5]),(Subscript[i, 3])]];
eqnu1=rearrangeComplexEqn[sqabs[(Subscript[u, 2])]==toMyForm[(Unom^2)]];
eqnu2=rearrangeComplexEqn[sqabs[(Subscript[u, 3])]==toMyForm[(Unom^2)]];
eqnSourceReactivePowerRatio=rearrangeComplexEqn[reactivePower[(Subscript[u, 1]),(Subscript[i, 1])]==1.5*reactivePower[(Subscript[u, 5]),(Subscript[i, 3])]];
eqnSettingAbsoluteTime={(Subscript[ure, 3])==0};
eqnsTogether=Union[eqnZ1,eqnZ2,eqnZ3,eqnZ4,eqnNod1,eqnNod2,eqnApparentPower1,eqnApparentPower2,eqnSourceActivePowerRatio1,eqnSourceActivePowerRatio2,eqnu1,eqnu2,eqnSourceReactivePowerRatio,eqnSettingAbsoluteTime];
unknownVariables=Union@Cases[eqnsTogether,Subscript[_,_],{0,∞}];
Length/@{eqnsTogether,unknownVariables}
vectfces=First/@eqnsTogether;
vars=unknownVariables;
grad[expr_]:=D[expr,#]&/@vars;
jacobian=grad/@vectfces;
mf@%;
xStart=unknownVariables/.{(Subscript[iim, index_]):>1,(Subscript[ire, index_]):>1,(Subscript[ure, index_]):>Unom,(Subscript[uim, index_]):>Unom};
xNext[xPrevious_]:=xPrevious-Inverse[jacobian/.Thread[vars->xPrevious]].(vectfces/.Thread[vars->xPrevious]);
aux=NestList[xNext,xStart,50][[-1]]

eqnVoltage3=rearrangeComplexEqn[sqabs[(Subscript[u, 3])]==toMyForm[(Unom^2)]];
eqnZ1=rearrangeComplexEqn[(Subscript[i, 1])==current[(Subscript[u, 1]),(Subscript[u, 3]),(Subscript[z, 1])]];
eqnZ2=rearrangeComplexEqn[(Subscript[i, 2])==current[(Subscript[u, 2]),(Subscript[u, 3]),(Subscript[z, 2])]];
eqnZ3=rearrangeComplexEqn[(Subscript[i, 3])==current[(Subscript[u, 3]),(Subscript[u, 4]),(Subscript[z, 3])]];
eqnNod3=rearrangeComplexEqn[(Subscript[i, 1])+(Subscript[i, 2])==(Subscript[i, 3])];
eqnNod4=rearrangeComplexEqn[(Subscript[i, 4])+(Subscript[i, 5])==(Subscript[i, 3])];
eqnSourceActivePowerRatio=rearrangeComplexEqn[activePower[(Subscript[u, 1]),(Subscript[i, 1])]==kReg*activePower[(Subscript[u, 2]),(Subscript[i, 2])]];
eqnSourcReactivePowerRatio=rearrangeComplexEqn[reactivePower[(Subscript[u, 1]),(Subscript[i, 1])]==kReg*reactivePower[(Subscript[u, 2]),(Subscript[i, 2])]];
eqnApparentPower4=rearrangeComplexEqn[apparentPower[(Subscript[u, 4]),(Subscript[i, 4])]==S4];
eqnApparentPower5=rearrangeComplexEqn[apparentPower[(Subscript[u, 4]),(Subscript[i, 5])]==S5];
eqnSettingAbsoluteTime={(Subscript[ure, 3])==0};
eqnsTogether=Union[eqnVoltage3,eqnZ1,eqnZ2,eqnZ3,eqnNod3,eqnNod4,eqnSourceActivePowerRatio,eqnSourcReactivePowerRatio,eqnApparentPower4,eqnApparentPower5,eqnSettingAbsoluteTime];
unknownVariables=Union@Cases[eqnsTogether,Subscript[_,_],{0,∞}];
Length/@{eqnsTogether,unknownVariables}

secondParamteterOfFindRoot=unknownVariables/.{(Subscript[iim, index_]):>{(Subscript[iim, index]),1},(Subscript[ire, index_]):>{(Subscript[ire, index]),1},(Subscript[ure, index_]):>{(Subscript[ure, index]),Unom},(Subscript[uim, index_]):>{(Subscript[uim, index]),Unom}};
solFindRoot=FindRoot[eqnsTogether,secondParamteterOfFindRoot]//Chop

mf=MatrixForm;
vectfces=First/@eqnsTogether;
vars=unknownVariables;
grad[expr_]:=D[expr,#]&/@vars;
jacobian=grad/@vectfces;
mf@%;
xStart=unknownVariables/.{(Subscript[iim, index_]):>1,(Subscript[ire, index_]):>1,(Subscript[ure, index_]):>Unom,(Subscript[uim, index_]):>Unom};
xNext[xPrevious_]:=xPrevious-Inverse[jacobian/.Thread[vars->xPrevious]].(vectfces/.Thread[vars->xPrevious]);
aux=NestList[xNext,xStart,50][[-1]];
Thread[vars->aux]
solFindRoot

mf@(Chop@jacobian)