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

(*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_]:=3times[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]}];

kz=100;
{(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]),(Subscript[z, 4]),(Subscript[z, 5]),(Subscript[z, 6])}=kz*(toMyForm/@{.01+0.05I,2*(01+0.05I),0.01+0.02I,0.01+0.02I,0.01+0.02I,0.01+0.02I});

powerFactor5=0.8;
S3=1*(10^6)*{1,Sin[ArcCos[1]]};
S5=5*(10^6)*{powerFactor5,Sin[ArcCos[powerFactor5]]};
{(Subscript[u, #])={(Subscript[ure, #]),(Subscript[uim, #])},(Subscript[i, #])={(Subscript[ire, #]),(Subscript[iim, #])}}&/@Range[10];
Unom=(22000/(Sqrt[3]));
kRegActive=1;
kRegReactive=2;

eqnVoltage5=rearrangeComplexEqn[sqabs[(Subscript[u, 5])]==toMyForm[(Unom^2)]];
eqnZ1=rearrangeComplexEqn[(Subscript[i, 1])==current[(Subscript[u, 1]),(Subscript[u, 4]),(Subscript[z, 1])]];
eqnZ2=rearrangeComplexEqn[(Subscript[i, 2])==current[(Subscript[u, 2]),(Subscript[u, 4]),(Subscript[z, 2])]];
eqnZ4=rearrangeComplexEqn[(Subscript[i, 5])==current[(Subscript[u, 3]),(Subscript[u, 6]),(Subscript[z, 4])]];
eqnZ3=rearrangeComplexEqn[(Subscript[i, 4])==current[(Subscript[u, 4]),(Subscript[u, 6]),(Subscript[z, 3])]];
eqnZ5=rearrangeComplexEqn[(Subscript[i, 6])==current[(Subscript[u, 6]),(Subscript[u, 5]),(Subscript[z, 5])]];
eqnZ6=rearrangeComplexEqn[(Subscript[i, 3])==current[(Subscript[u, 4]),(Subscript[u, 5]),(Subscript[z, 6])]];
eqnNod4=rearrangeComplexEqn[(Subscript[i, 1])+(Subscript[i, 2])==(Subscript[i, 3])+(Subscript[i, 4])];
eqnNod6=rearrangeComplexEqn[(Subscript[i, 5])+(Subscript[i, 4])==(Subscript[i, 6])];
eqnNod5=rearrangeComplexEqn[(Subscript[i, 6])+(Subscript[i, 3])==(Subscript[i, 7])];
eqnApparentPower5=rearrangeComplexEqn[apparentPower[(Subscript[u, 5]),(Subscript[i, 7])]==S5];
eqnApparentPower3=rearrangeComplexEqn[apparentPower[(Subscript[u, 3]),(Subscript[i, 5])]==S5];
eqnSettingAbsoluteTime={(Subscript[ure, 6])==0};
eqnSourceActivePowerRatio=rearrangeComplexEqn[activePower[(Subscript[u, 1]),(Subscript[i, 1])]==kRegActive*activePower[(Subscript[u, 2]),(Subscript[i, 2])]];
eqnSourcReactivePowerRatio=rearrangeComplexEqn[reactivePower[(Subscript[u, 1]),(Subscript[i, 1])]==kRegReactive*reactivePower[(Subscript[u, 2]),(Subscript[i, 2])]];
eqnsTogether=Union[eqnVoltage5,eqnZ1,eqnZ2,eqnZ4,eqnZ3,eqnZ5,eqnZ6,eqnNod4,eqnNod6,eqnNod5,eqnApparentPower5,eqnSettingAbsoluteTime,eqnApparentPower3,eqnSourceActivePowerRatio,eqnSourcReactivePowerRatio];
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

Inverse[(Chop@jacobian)]