diff --git a/TODO.md b/TODO.md
index f39754f2c..89d831bef 100644
--- a/TODO.md
+++ b/TODO.md
@@ -10,8 +10,8 @@
 - [x] Improve the speed of `DAETimeSeries.unpack`.
 - [x] Allow simulating to a time, pause, and continue to a new ending time.
 - [x] Store TDS data as NumPy compressed format `npz`; allow to reload both `npy` and `npz`.
+- [x] Power flow model variables erroneously point to the old `dae.y` array (fixed)
 *   Allow loading default config for `selftest`
-*   Power flow model variables erroneously point to the old `dae.y` array
 
 ## Later Versions
 *   Generalize two-terminal and multi-terminal devices
diff --git a/andes/models/line.py b/andes/models/line.py
index f979ce461..3189d698d 100644
--- a/andes/models/line.py
+++ b/andes/models/line.py
@@ -44,6 +44,7 @@ def __init__(self, system=None, config=None):
         Model.__init__(self, system, config)
         self.group = 'ACLine'
         self.flags.pflow = True
+        self.flags.tds = True
 
         self.a1 = ExtAlgeb(model='Bus', src='a', indexer=self.bus1, tex_name='a_1',
                            info='phase angle of the from bus')
diff --git a/andes/routines/tds.py b/andes/routines/tds.py
index 04e51fa63..a73de7b4f 100644
--- a/andes/routines/tds.py
+++ b/andes/routines/tds.py
@@ -61,7 +61,6 @@ def __init__(self, system=None, config=None):
         self.deltatmax = 0
         self.h = 0
         self.next_pc = 0
-        self.eye = None
         self.Teye = None
         self.qg = np.array([])
         self.tol_zero = self.config.tol / 100
@@ -100,17 +99,25 @@ def init(self):
 
         self._reset()
         self._load_pert()
-        system.set_address(models=system.exist.tds)
+
+        # Note:
+        #   calling `set_address` on `system.exist.pflow_tds` will point all variables
+        #   to the new array after extending `dae.y`
+        system.set_address(models=system.exist.pflow_tds)
+
         system.set_dae_names(models=system.exist.tds)
 
         system.dae.clear_ts()
         system.store_sparse_pattern(models=system.exist.pflow_tds)
         system.store_adder_setter(models=system.exist.pflow_tds)
         system.vars_to_models()
+
+        # Initialize `system.exist.tds` only to avoid Bus overwriting power flow solutions
         system.init(system.exist.tds)
         system.store_switch_times(system.exist.tds)
-        self.eye = spdiag([1] * system.dae.n)
-        self.Teye = spdiag(system.dae.Tf.tolist()) * self.eye
+
+        # Build mass matrix into `self.Teye`
+        self.Teye = spdiag(system.dae.Tf.tolist())
         self.qg = np.zeros(system.dae.n + system.dae.m)
 
         self.initialized = self.test_init()
@@ -708,7 +715,6 @@ def _reset(self):
         self.deltatmax = 0
         self.h = 0
         self.next_pc = 0.1
-        self.eye = None
         self.Teye = None
         self.qg = np.array([])