diff --git a/docs/whats_new/v0-9-14.rst b/docs/whats_new/v0-9-14.rst
index 45a1afd4c..bd641d520 100644
--- a/docs/whats_new/v0-9-14.rst
+++ b/docs/whats_new/v0-9-14.rst
@@ -14,7 +14,10 @@ Bug Fixes
   silent error because if not used, it did not do any trouble. Only when
   importing a :code:`Network` and reconstructing the component parameters the
   :code:`CharMap` could not be constructed successfully because it was missing
-  dimensions (`PR #934 <https://github.com/oemof/tespy/pull/934>`__).
+  dimensions (`PR #933 <https://github.com/oemof/tespy/pull/933>`__).
+- The reactor classes :code:`FuelCell` and :code:`WaterElectrolyzer` were
+  missing the implementation for a :code:`PowerConnection`
+  (`PR #934 <https://github.com/oemof/tespy/pull/934>`__).
 
 Other Changes
 #############
diff --git a/src/tespy/components/reactors/fuel_cell.py b/src/tespy/components/reactors/fuel_cell.py
index 9f7e8ff05..c9f51e49c 100644
--- a/src/tespy/components/reactors/fuel_cell.py
+++ b/src/tespy/components/reactors/fuel_cell.py
@@ -214,7 +214,7 @@ def get_parameters(self):
         }
 
     def get_mandatory_constraints(self):
-        return {
+        constraints = {
             'mass_flow_constraints': dc_cmc(**{
                 'func': self.reactor_mass_flow_func,
                 'deriv': self.reactor_mass_flow_deriv,
@@ -247,6 +247,14 @@ def get_mandatory_constraints(self):
                 "description": "reactor pressure equality equations"
             })
         }
+        if len(self.power_outl) > 0:
+            constraints["energy_connector_balance"] = dc_cmc(**{
+                "func": self.energy_connector_balance_func,
+                "dependents": self.energy_connector_dependents,
+                "num_eq_sets": 1
+            })
+
+        return constraints
 
     @staticmethod
     def get_bypass_constraints():
@@ -260,6 +268,10 @@ def inlets():
     def outlets():
         return ['out1', 'out2']
 
+    @staticmethod
+    def poweroutlets():
+        return ["power"]
+
     def _add_missing_fluids(self, connections):
         if self.inl[1] in connections:
             return ["O2"]
@@ -325,6 +337,28 @@ def calc_e0(self):
 
         return e0
 
+    def energy_connector_balance_func(self):
+        r"""
+        (optional) energy balance equation connecting the power connector to
+        the component's power. Since the power of the FuelCell is negative by
+        definition of the system, the two quantities are added in the residual
+        so the power on the PowerConnection is positive again.
+
+        Returns
+        -------
+        residual : float
+            Residual value of equation
+
+            .. math::
+
+                0=\dot E + P
+        """
+        return self.power_outl[0].E.val_SI + self.P.val_SI
+
+    def energy_connector_dependents(self):
+        return [self.power_outl[0].E, self.P]
+
+
     def eta_func(self):
         r"""
         Equation for efficiency.
@@ -454,7 +488,7 @@ def energy_balance_deriv(self, increment_filter, k, dependents=None):
 
         # derivatives for variable P
         if self.P.is_var:
-            self.jacobian[k, self.p.J_col] = 1
+            self.jacobian[k, self.P.J_col] = 1
 
     def energy_balance_dependents(self):
         return [
diff --git a/src/tespy/components/reactors/water_electrolyzer.py b/src/tespy/components/reactors/water_electrolyzer.py
index 8235a0c6c..c7da6b4a4 100644
--- a/src/tespy/components/reactors/water_electrolyzer.py
+++ b/src/tespy/components/reactors/water_electrolyzer.py
@@ -253,7 +253,7 @@ def get_parameters(self):
         }
 
     def get_mandatory_constraints(self):
-        return {
+        constraints = {
             'mass_flow_constraints': dc_cmc(**{
                 'func': self.reactor_mass_flow_func,
                 'deriv': self.reactor_mass_flow_deriv,
@@ -291,6 +291,14 @@ def get_mandatory_constraints(self):
                 "description": "equation for same temperature of product gases"
             })
         }
+        if len(self.power_inl) > 0:
+            constraints["energy_connector_balance"] = dc_cmc(**{
+                "func": self.energy_connector_balance_func,
+                "dependents": self.energy_connector_dependents,
+                "num_eq_sets": 1
+            })
+
+        return constraints
 
     @staticmethod
     def get_bypass_constraints():
@@ -304,6 +312,10 @@ def inlets():
     def outlets():
         return ['out1', 'out2', 'out3']
 
+    @staticmethod
+    def powerinlets():
+        return ["power"]
+
     def _add_missing_fluids(self, connections):
         if self.inl[1] in connections:
             return ["H2O"]
@@ -369,6 +381,25 @@ def calc_e0(self):
 
         return e0
 
+    def energy_connector_balance_func(self):
+        r"""
+        (optional) energy balance equation connecting the power connector to
+        the component's power
+
+        Returns
+        -------
+        residual : float
+            Residual value of equation
+
+            .. math::
+
+                0=\dot E - P
+        """
+        return self.power_inl[0].E.val_SI - self.P.val_SI
+
+    def energy_connector_dependents(self):
+        return [self.power_inl[0].E, self.P]
+
     def gas_temperature_func(self):
         r"""
         Equation for temperature equality of product gases.
diff --git a/tests/test_components/test_reactors.py b/tests/test_components/test_reactors.py
index 7f5ecff92..36b74922f 100644
--- a/tests/test_components/test_reactors.py
+++ b/tests/test_components/test_reactors.py
@@ -9,16 +9,21 @@
 
 SPDX-License-Identifier: MIT
 """
+from pytest import approx
 
+from tespy.components import FuelCell
+from tespy.components import PowerSink
+from tespy.components import PowerSource
 from tespy.components import Sink
 from tespy.components import Source
 from tespy.components import WaterElectrolyzer
 from tespy.connections import Bus
 from tespy.connections import Connection
+from tespy.connections import PowerConnection
 from tespy.networks import Network
 
 
-class TestReactors:
+class TestWaterElectrolyzer:
 
     def setup_method(self):
         """Set up network for electrolyzer tests."""
@@ -64,16 +69,16 @@ def test_WaterElectrolyzer(self, tmp_path):
         self.nw.solve('design')
         self.nw.assert_convergence()
         assert self.nw.status == 0
-        msg = ('Value of power must be ' + str(power.P.val) + ', is ' +
-               str(self.instance.P.val) + '.')
-        assert round(power.P.val, 1) == round(self.instance.P.val), msg
+        msg = (
+            f"Value of power must be {power.P.val}, is {self.instance.P.val}."
+        )
+        assert approx(power.P.val) == self.instance.P.val, msg
 
         # effieciency was set to 100 % with inlet and outlet states of the
         # reaction educts and products beeing identical to reference state
         # therefore Q must be equal to 0
-        msg = ('Value of heat output must be 0.0, is ' +
-               str(self.instance.Q.val) + '.')
-        assert round(self.instance.Q.val, 4) == 0.0, msg
+        msg = f"Value of heat must be 0.0, is {self.instance.Q.val}."
+        assert approx(self.instance.Q.val, abs=1e-4) == 0.0, msg
 
         # reset power, change efficiency value and specify heat bus value
         power.set_attr(P=None)
@@ -89,9 +94,10 @@ def test_WaterElectrolyzer(self, tmp_path):
 
         self.nw.solve('design')
         self.nw.assert_convergence()
-        msg = ('Value of heat flow must be ' + str(heat.P.val) +
-               ', is ' + str(self.instance.Q.val) + '.')
-        assert round(heat.P.val, 1) == round(self.instance.Q.val), msg
+        msg = (
+            f"Value of heat must be {heat.P.val}, is {self.instance.Q.val}."
+        )
+        assert approx(heat.P.val) == self.instance.Q.val, msg
         self.nw.save(tmp_path)
 
         # check bus function:
@@ -100,9 +106,8 @@ def test_WaterElectrolyzer(self, tmp_path):
         heat.set_attr(P=Q)
         self.nw.solve('offdesign', design_path=tmp_path)
         self.nw.assert_convergence()
-        msg = ('Value of heat flow must be ' + str(Q) +
-               ', is ' + str(self.instance.Q.val) + '.')
-        assert round(Q, 1) == round(self.instance.Q.val), msg
+        msg = f"Value of heat must be {Q}, is {self.instance.Q.val}."
+        assert approx(Q) == self.instance.Q.val, msg
 
         # delete both busses again
         self.nw.del_busses(heat, power)
@@ -112,11 +117,11 @@ def test_WaterElectrolyzer(self, tmp_path):
         self.instance.set_attr(eta=0.9, e='var')
         self.nw.solve('design')
         self.nw.assert_convergence()
-        msg = ('Value of efficiency must be ' + str(self.instance.eta.val) +
-               ', is ' + str(self.instance.e0 / self.instance.e.val) + '.')
-        eta = round(self.instance.eta.val, 2)
-        eta_calc = round(self.instance.e0 / self.instance.e.val, 2)
-        assert eta == eta_calc, msg
+        msg = (
+            f"Value of efficiency must be {self.instance.eta.val}, is "
+            f"{self.instance.e0 / self.instance.e.val}."
+        )
+        assert approx(self.instance.eta.val) == self.instance.e0 / self.instance.e.val, msg
 
         # test efficiency value > 1, Q must be larger than 0
         e = 130e6
@@ -125,14 +130,16 @@ def test_WaterElectrolyzer(self, tmp_path):
         self.nw.solve('design')
         self.nw.assert_convergence()
         # test efficiency
-        msg = ('Value of efficiency must be ' + str(self.instance.e0 / e) +
-               ', is ' + str(self.instance.eta.val) + '.')
-        eta = round(self.instance.e0 / e, 2)
-        eta_calc = round(self.instance.eta.val, 2)
-        assert eta == eta_calc, msg
+        msg = (
+            f"Value of efficiency must be {self.instance.e0 / e}, is "
+            f"{self.instance.eta.val}."
+        )
+        assert approx(self.instance.e0 / e) == self.instance.eta.val, msg
         # test Q
-        msg = ('Value of heat must be larger than zero, is ' +
-               str(self.instance.Q.val) + '.')
+        msg = (
+            f"Value of heat must be larger than zero, is "
+            f"{self.instance.Q.val}."
+        )
         assert self.instance.Q.val > 0, msg
 
         # test specific energy consumption
@@ -141,9 +148,11 @@ def test_WaterElectrolyzer(self, tmp_path):
         self.instance.set_attr(e=e)
         self.nw.solve('design')
         self.nw.assert_convergence()
-        msg = ('Value of specific energy consumption e must be ' + str(e) +
-               ', is ' + str(self.instance.e.val) + '.')
-        assert round(e, 1) == round(self.instance.e.val, 1), msg
+        msg = (
+            f"Value of specific energy consumption e must be {e}, is "
+            f"{self.instance.e.val}."
+        )
+        assert approx(e) == self.instance.e.val, msg
 
         # test cooling loop pressure ratio, zeta as variable value
         pr = 0.95
@@ -152,24 +161,156 @@ def test_WaterElectrolyzer(self, tmp_path):
         self.nw.solve('design')
         self.nw.save(tmp_path)
         self.nw.assert_convergence()
-        msg = ('Value of pressure ratio must be ' + str(pr) + ', is ' +
-               str(self.instance.pr.val) + '.')
-        assert round(pr, 2) == round(self.instance.pr.val, 2), msg
+        msg = (
+            f"Value of pressure ratio must be {pr}, is {self.instance.pr.val}."
+        )
+        assert approx(pr) == self.instance.pr.val, msg
 
         # use zeta as offdesign parameter, at design point pressure
         # ratio must not change
         self.instance.set_attr(zeta=None, offdesign=['zeta'])
         self.nw.solve('offdesign', design_path=tmp_path)
         self.nw.assert_convergence()
-        msg = ('Value of pressure ratio must be ' + str(pr) + ', is ' +
-               str(self.instance.pr.val) + '.')
-        assert round(pr, 2) == round(self.instance.pr.val, 2), msg
+        msg = (
+            f"Value of pressure ratio must be {pr}, is {self.instance.pr.val}."
+        )
+        assert approx(pr) == self.instance.pr.val, msg
 
         # test heat output specification in offdesign mode
         Q = self.instance.Q.val * 0.9
         self.instance.set_attr(Q=Q, P=None)
         self.nw.solve('offdesign', design_path=tmp_path)
         self.nw.assert_convergence()
-        msg = ('Value of heat must be ' + str(Q) + ', is ' +
-               str(self.instance.Q.val) + '.')
-        assert round(Q, 0) == round(self.instance.Q.val, 0), msg
+        msg = f"Value of heat must be {Q}, is {self.instance.Q.val}."
+        assert approx(Q) == self.instance.Q.val, msg
+
+    def test_WaterElectrolyzer_PowerConnection(self):
+        """test utilization of PowerConnection with WaterElectrolyzer"""
+        self.nw.get_conn('h2o').set_attr(T=25, p=1)
+        self.nw.get_conn('h2').set_attr(T=25)
+
+        power = PowerSource("grid")
+        e1 = PowerConnection(power, "power", self.instance, "power", label="e1")
+        self.nw.add_conns(e1)
+
+        e1.set_attr(E=2.5e6)
+        self.nw.solve('design')
+
+        assert approx(self.instance.P.val_SI) == e1.E.val_SI
+
+
+class TestFuelCell:
+
+    def setup_method(self):
+        """Set up network for fuel cell tests."""
+        self.nw = Network()
+        self.nw.units.set_defaults(**{
+            "pressure": "bar", "temperature": "degC"
+        })
+        self.instance = FuelCell('fuel cell')
+
+        cw_in = Source('cooling water in')
+        cw_out = Sink('cooling water out')
+        o2_source = Source('oxygen source')
+        h2_source = Source('hydrogen source')
+        water_sink = Sink('water sink')
+
+        self.instance.set_attr()
+
+        cw_fc = Connection(
+            cw_in, 'out1', self.instance, 'in1', fluid={'H2O': 1}, T=20, p=1
+        )
+        fc_cw = Connection(self.instance, 'out1', cw_out, 'in1', T=40)
+        self.nw.add_conns(cw_fc, fc_cw)
+
+        o2_fc = Connection(
+            o2_source, 'out1', self.instance, 'in2', label='o2', T=25, p=1
+        )
+        h2_fc = Connection(
+            h2_source, 'out1', self.instance, 'in3', label='h2', T=25
+        )
+        fc_h2o = Connection(
+            self.instance, 'out2', water_sink, 'in1', label='h2o', T=50
+        )
+        self.nw.add_conns(o2_fc, h2_fc, fc_h2o)
+
+    def test_FuelCell(self, tmp_path):
+        """Test component properties of fuel cell."""
+        tmp_path = f'{tmp_path}.json'
+
+        # test efficiency vs. specific energy consumption
+        # eta = e / e0, with both e and e0 negative for the fuel cell
+        self.nw.get_conn('h2').set_attr(m=0.01)
+        self.instance.set_attr(pr=0.98)
+        self.instance.set_attr(eta=0.45)
+        self.nw.solve('design')
+        self.nw.assert_convergence()
+        eta = self.instance.eta.val
+        eta_calc = self.instance.e.val / self.instance.e0
+        msg = f"Value of efficiency must be {eta}, is {eta_calc}."
+        assert approx(eta) == eta_calc, msg
+
+        # at eta < 1 the fuel cell rejects heat to the cooling loop (Q < 0)
+        msg = (
+            f"Value of heat must be less than zero, is "
+            f"{self.instance.Q.val}."
+        )
+        assert self.instance.Q.val < 0, msg
+
+        # test specific energy consumption specification
+        e = self.instance.e0 * 0.5
+        self.instance.set_attr(e=None, eta=None)
+        self.instance.set_attr(e=e)
+        self.nw.solve('design')
+        self.nw.assert_convergence()
+        msg = (
+            f"Value of specific energy consumption e must be {e}, is "
+            f"{self.instance.e.val}."
+        )
+        assert approx(e) == self.instance.e.val, msg
+
+        # test cooling loop pressure ratio, zeta as variable value
+        pr = 0.95
+        self.instance.set_attr(
+            pr=pr, e=None, eta=None, zeta='var', P=-2e5, design=['pr']
+        )
+        self.nw.solve('design')
+        self.nw.save(tmp_path)
+        self.nw.assert_convergence()
+        msg = (
+            f"Value of pressure ratio must be {pr}, is {self.instance.pr.val}."
+        )
+        assert approx(pr) == self.instance.pr.val, msg
+
+        # use zeta as offdesign parameter, at design point pressure
+        # ratio must not change
+        self.instance.set_attr(zeta=None, offdesign=['zeta'])
+        self.nw.solve('offdesign', design_path=tmp_path)
+        self.nw.assert_convergence()
+        msg = (
+            f"Value of pressure ratio must be {pr}, is {self.instance.pr.val}."
+        )
+        assert approx(pr) == self.instance.pr.val, msg
+
+        # test heat output specification in offdesign mode
+        Q = self.instance.Q.val * 0.9
+        self.instance.set_attr(Q=Q, P=None)
+        self.nw.solve('offdesign', design_path=tmp_path)
+        self.nw.assert_convergence()
+        msg = f"Value of heat must be {Q}, is {self.instance.Q.val}."
+        assert approx(Q) == self.instance.Q.val, msg
+
+    def test_FuelCell_PowerConnection(self):
+        """Test utilization of PowerConnection with FuelCell."""
+        self.instance.set_attr(eta=.5, pr=0.98)
+        power_sink = PowerSink("grid")
+        e1 = PowerConnection(
+            self.instance, "power", power_sink, "power", label="e1"
+        )
+        self.nw.add_conns(e1)
+
+        e1.set_attr(E=2e5)
+        self.nw.solve('design')
+        self.nw.assert_convergence()
+
+        assert approx(-self.instance.P.val_SI) == e1.E.val_SI