Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];

c=4186.;
mHW=0.5;
m2=2.;
m3=1;
Ploses=1000;
COP=0.5*((T4+273)/(T4-T7));
T1=5;
T4=55;
T6=35;

kS=5000;
(*heat exchanger*)
Δ1=T6-T2;
Δ2=T7-T1;
eqnsExchanger={m2*c*T6+mHW*c*T1==m2*c*T7+mHW*c*T2,
QkS==kS*((Δ1-Δ2)/(Log[((Δ1)/(Δ2))])),
QkS==m2*c*(T6-T7)
};
eqnsExchangerLinear={m2*c*T6+mHW*c*T1==m2*c*T7+mHW*c*T2,
QkS==kS*((Δ1+Δ2)/2),
QkS==m2*c*(T6-T7)
};
(*heat pump*)
eqnsHeatPump={m2*c*T7+m3*c*T3+Pel==m2*c*T8+m3*c*T4,
QHP==m3*c*(T4-T3),
QHP==Pel*COP
};
(*heat loses*)
eqnHeatLoses={m3*c*T4==Ploses+m3*c*T5};
(*mixing*)
eqnsMixing={(m3-mHW)*c*T5+mHW*c*T2==m3*c*T3};
eqnsTogether=Union[eqnsExchanger,eqnsHeatPump,eqnsMixing,eqnHeatLoses];
Union[Cases[eqnsTogether,_Symbol,All]]//Length;
eqnsTogether//Length;
eqnsTogetherLinear=Union[eqnsExchangerLinear,eqnsHeatPump,eqnsMixing,eqnHeatLoses];
solLinear=Quiet@Solve[eqnsTogetherLinear][[1]];
start=solLinear/.Rule->List;
solFinal=FindRoot[eqnsTogether,start]

powerObtained=QHP+QkS/.solFinal;
powerElectro=Pel/.solFinal;
Print["we obtained ",(10^(-3))powerObtained," kW  and we pay for ",(10^(-3))powerElectro," kW"]