dongxuelian2 · Oct 7, 2026
diff --git a/control/stochsys.py b/control/stochsys.py
@@ -83,6 +83,10 @@ def lqe(*args, symmetric_kwargs=None, **kwargs):
         Set the method used for computing the result.  Current methods are
         'slycot' and 'scipy'.  If set to None (default), try 'slycot' first
         and then 'scipy'.
+    return_filter_form : bool, optional
+        For a discrete-time `sys`, return the measurement-update gain instead
+        of the predictor gain.  Default is False.  See `dlqe` for details.
+        Only supported for discrete-time systems.
     symmetric_kwargs : dict, optional
         Keyword arguments passed to `scipy.linalg.issymmetric` or
         `scipy.linalg.ishermitian`.
@@ -132,7 +136,7 @@ def lqe(*args, symmetric_kwargs=None, **kwargs):
     # If we were passed a discrete-time system as the first arg, use dlqe()
     if isinstance(args[0], LTI) and isdtime(args[0], strict=True):
         # Call dlqe
-        return dlqe(*args, **kwargs)
+        return dlqe(*args, symmetric_kwargs=symmetric_kwargs, **kwargs)
 
     # Get the method to use (if specified as a keyword)
     method = kwargs.pop('method', None)
@@ -186,8 +190,8 @@ def lqe(*args, symmetric_kwargs=None, **kwargs):
 
 
 # contributed by Sawyer B. Fuller <minster@uw.edu>
-def dlqe(*args, symmetric_kwargs=None, **kwargs):
-    r"""dlqe(A, G, C, QN, RN, [, N])
+def dlqe(*args, return_filter_form=False, symmetric_kwargs=None, **kwargs):
+    r"""dlqe(A, G, C, QN, RN, [, NN])
 
     Discrete-time linear quadratic estimator (Kalman filter).
 
@@ -207,14 +211,31 @@ def dlqe(*args, symmetric_kwargs=None, **kwargs):
 
     .. math:: x_e[n+1] = A x_e[n] + B u[n] + L(y[n] - C x_e[n] - D u[n])
 
-    produces a state estimate x_e[n] that minimizes the mean squared
-    estimation error x[n] - x_e[n] using the sensor measurements y. The
-    noise cross-correlation `NN` is set to zero when omitted.
+    produces a prior state estimate x_e[n] using measurements through y[n-1].
+    This estimate minimizes the mean squared estimation error x[n] - x_e[n].
+    The noise cross-correlation `NN` is set to zero when omitted.
+
+    If `return_filter_form` is True, the returned gain is instead the
+    measurement-update gain L_f for the current state estimate:
+
+    .. math::
+
+        x_f[n] &= x_e[n] + L_f (y[n] - C x_e[n] - D u[n]) \\
+        L_f &= P C^T (C P C^T + RN)^{-1}
+
+    The next prior estimate is A x_f[n] + B u[n].  The predictor gain is
+    L = A L_f, but computing L_f does not require A to be invertible.
+
+    The system matrices can also be supplied as a discrete-time `StateSpace`
+    system: ``L, P, E = dlqe(sys, QN, RN)``.
 
     Parameters
     ----------
     A, G, C : 2D array_like
         Dynamics, process noise (disturbance), and output matrices.
+    sys : `StateSpace`
+        Discrete-time linear I/O system, with the process noise input taken
+        as the system input.
     QN, RN : 2D array_like
         Process and sensor noise covariance matrices.
     NN : 2D array, optional
@@ -223,28 +244,40 @@ def dlqe(*args, symmetric_kwargs=None, **kwargs):
         Set the method used for computing the result.  Current methods are
         'slycot' and 'scipy'.  If set to None (default), try 'slycot'
         first and then 'scipy'.
+    return_filter_form : bool, optional
+        Return the measurement-update gain L_f if True.  If False (default),
+        return the predictor gain L = A L_f.  The covariance `P` and poles `E`
+        are the same for both forms.
     symmetric_kwargs : dict, optional
         Keyword arguments passed to `scipy.linalg.issymmetric` or
         `scipy.linalg.ishermitian`.
 
     Returns
     -------
     L : 2D array
-        Kalman estimator gain.
+        Kalman predictor gain, or measurement-update gain if
+        `return_filter_form` is True.
     P : 2D array
-        Solution to Riccati equation.
+        Prior estimation error covariance (before the measurement update),
+        given by the solution to the discrete-time Riccati equation:
 
         .. math::
 
-            A P + P A^T - (P C^T + G N) R^{-1}  (C P + N^T G^T) + G Q G^T = 0
+            P &= A P A^T + G QN G^T \\
+              &\quad - A P C^T (C P C^T + RN)^{-1} C P A^T
 
     E : 1D array
-        Eigenvalues of estimator poles eig(A - L C).
+        Eigenvalues of the estimator error dynamics: eig(A - L C) for the
+        predictor form, or eig(A - A L C) for the filter form.
 
     Examples
     --------
     >>> L, P, E = dlqe(A, G, C, QN, RN)                         # doctest: +SKIP
     >>> L, P, E = dlqe(A, G, C, QN, RN, NN)                     # doctest: +SKIP
+    >>> import control as ct
+    >>> M, P, E = ct.dlqe(0, 1, 1, 1, 1, return_filter_form=True)
+    >>> M
+    array([[0.5]])
 
     See Also
     --------
@@ -306,6 +339,10 @@ def dlqe(*args, symmetric_kwargs=None, **kwargs):
     # Compute the result (dimension and symmetry checking done in dare())
     P, E, LT = dare(A.T, C.T, G @ QN @ G.T, RN, method=method,
                     _Bs="C", _Qs="QN", _Rs="RN", _Ss="NN", symmetric_kwargs=symmetric_kwargs)
+    if return_filter_form:
+        # Solve for the measurement-update gain using the prior covariance.
+        # This also works when A is singular.
+        return np.linalg.solve(C @ P @ C.T + RN, C @ P).T, P, E
     return LT.T, P, E
 
 

diff --git a/control/tests/stochsys_test.py b/control/tests/stochsys_test.py
@@ -59,7 +59,9 @@ def test_lqe_symmetric_kwargs():
         symmetric_kwargs={"rtol": 1e-12},
     )
 
-def test_dlqe_symmetric_kwargs():
+@pytest.mark.parametrize("return_filter_form", [False, True])
+@pytest.mark.parametrize("cdlqe", [lqe, dlqe])
+def test_dlqe_symmetric_kwargs(cdlqe, return_filter_form):
     A = np.array([
         [0.5, 0.],
         [0., 0.4]
@@ -71,14 +73,16 @@ def test_dlqe_symmetric_kwargs():
         [0.2 + 1e-15, 1.]
     ])
     RN = np.eye(2)
+    sys = ct.ss(A, G, C, 0, dt=True)
 
     # Exact symmetry check should fail
     with pytest.raises(ControlArgument, match="symmetric"):
-        dlqe(A, G, C, QN, RN, method="scipy")
+        cdlqe(sys, QN, RN, method="scipy",
+              return_filter_form=return_filter_form)
 
     # Passing rtol should allow the nearly symmetric matrix
-    dlqe(
-        A, G, C, QN, RN,
+    cdlqe(
+        sys, QN, RN, return_filter_form=return_filter_form,
         method="scipy",
         symmetric_kwargs={"rtol": 1e-12},
     )
@@ -133,6 +137,100 @@ def test_DLQE(method):
     L, P, poles = dlqe(A, G, C, QN, RN, method=method)
     check_DLQE(L, P, poles, G, QN, RN)
 
+
+@pytest.mark.parametrize("method", [
+    None, pytest.param('slycot', marks=pytest.mark.slycot), 'scipy'])
+@pytest.mark.parametrize("a", [0., 0.7, 1.2])
+@pytest.mark.parametrize("return_filter_form", [False, True])
+def test_dlqe_scalar_gain(method, a, return_filter_form):
+    # The scalar DARE is P**2 + (r * (1 - a**2) - q) * P - q * r = 0.
+    g, qn, r = 0.3, 2., 0.7
+    q = g**2 * qn
+    b = r * (1 - a**2) - q
+    p = (-b + np.sqrt(b**2 + 4 * q * r)) / 2
+    m = p / (p + r)
+
+    L, P, E = dlqe(a, g, 1, qn, r, method=method,
+                   return_filter_form=return_filter_form)
+    assert L.shape == P.shape == (1, 1)
+    assert E.shape == (1,)
+    np.testing.assert_allclose(P, [[p]], rtol=1e-10)
+    np.testing.assert_allclose(L, [[m if return_filter_form else a * m]],
+                               rtol=1e-10, atol=1e-14)
+    np.testing.assert_allclose(E, [a * (1 - m)], rtol=1e-10, atol=1e-14)
+
+
+@pytest.mark.parametrize("method", [
+    None, pytest.param('slycot', marks=pytest.mark.slycot), 'scipy'])
+@pytest.mark.parametrize("A", [
+    [[1.1, 0.2, 0], [0, 0.7, 0.1], [0, 0, 0.3]],
+    [[1.1, 0.2, 0], [0, 0.7, 0.1], [0, 0, 0]],
+])
+def test_dlqe_filter_form_covariance(A, method):
+    A = np.array(A)
+    G = np.array([[1, 0.25], [0.5, 1], [1, -0.5]])
+    C = np.array([[1, 0, 0.25], [0.5, 1, 0]])
+    QN = np.array([[0.5, 0.125], [0.125, 0.75]])
+    RN = np.array([[0.75, 0.25], [0.25, 1.5]])
+
+    # Independently converge the measurement update and prediction, using
+    # the Joseph form to compute the posterior covariance.
+    prior = np.eye(3)
+    for _ in range(200):
+        M = np.linalg.solve(C @ prior @ C.T + RN, C @ prior).T
+        correction = np.eye(3) - M @ C
+        posterior = correction @ prior @ correction.T + M @ RN @ M.T
+        next_prior = A @ posterior @ A.T + G @ QN @ G.T
+        if np.linalg.norm(next_prior - prior) < 1e-13:
+            break
+        prior = next_prior
+    else:
+        pytest.fail("Kalman covariance iteration did not converge")
+
+    L, P, E = dlqe(A, G, C, QN, RN, method=method)
+    Lpred, Ppred, Epred = dlqe(A, G, C, QN, RN, method=method,
+                              return_filter_form=False)
+    Lfilter, Pfilter, Efilter = dlqe(A, G, C, QN, RN, method=method,
+                                    return_filter_form=True)
+    np.testing.assert_allclose(P, prior, rtol=1e-10, atol=1e-12)
+    np.testing.assert_allclose(Lfilter, M, rtol=1e-10, atol=1e-12)
+    np.testing.assert_allclose(L, A @ M, rtol=1e-10, atol=1e-12)
+    np.testing.assert_allclose(Lpred, L)
+    np.testing.assert_allclose(Ppred, P)
+    np.testing.assert_allclose(Pfilter, P)
+    np.testing.assert_allclose(Epred, E)
+    np.testing.assert_allclose(Efilter, E)
+    np.testing.assert_allclose(
+        np.sort_complex(E), np.sort_complex(np.linalg.eigvals(A @ correction)),
+        rtol=1e-10, atol=1e-12)
+    assert np.max(np.abs(E)) < 1
+
+    # The two-step filter and the existing predictor give the same next
+    # estimate, even when A is singular and cannot recover M from A @ M.
+    xhat, innovation = np.array([0.1, -0.2, 0.3]), np.array([0.4, -0.5])
+    np.testing.assert_allclose(
+        A @ (xhat + Lfilter @ innovation), A @ xhat + L @ innovation)
+    np.testing.assert_allclose(
+        P, A @ posterior @ A.T + G @ QN @ G.T, rtol=1e-10, atol=1e-12)
+
+
+@pytest.mark.parametrize("return_filter_form", [False, True])
+@pytest.mark.parametrize("dt", [True, 0.1])
+def test_dlqe_filter_form_call_format(return_filter_form, dt):
+    sys = ct.ss([[0.7, 0.2], [0, 0]], [[1], [0.3]], [[1, 0.4]], 0, dt=dt)
+    expected = dlqe(sys.A, sys.B, sys.C, 0.5, 0.8,
+                    return_filter_form=return_filter_form)
+    for cdlqe in (dlqe, lqe):
+        result = cdlqe(sys, 0.5, 0.8, return_filter_form=return_filter_form)
+        for actual, reference in zip(result, expected):
+            np.testing.assert_allclose(actual, reference)
+        with pytest.raises(ct.ControlNotImplemented, match="cross-covariance"):
+            cdlqe(sys, 0.5, 0.8, 0, return_filter_form=return_filter_form)
+        with pytest.raises(TypeError, match="unrecognized keyword"):
+            cdlqe(sys, 0.5, 0.8, return_filter_form=return_filter_form,
+                  unknown=True)
+
+
 def test_lqe_discrete():
     """Test overloading of lqe operator for discrete-time systems"""
     csys = ct.rss(2, 1, 1)

diff --git a/doc/stochastic.rst b/doc/stochastic.rst
@@ -182,10 +182,42 @@ called in several forms:
 
 where :code:`sys` is an :class:`LTI` object, and `A`, `G`, `C`, `QN`, `RN`,
 and `NN` are 2D arrays of appropriate dimension.  If :code:`sys` is a
-discrete-time system, the first two forms will compute the discrete
-time optimal controller.  For the second two forms, the :func:`dlqr`
+discrete-time system, the first two forms will compute the discrete-time
+optimal estimator.  For the second two forms, the :func:`dlqe`
 function can be used.  Additional arguments and details are given on
-the :func:`lqr` and :func:`dlqr` documentation pages.
+the :func:`lqe` and :func:`dlqe` documentation pages.
+
+In discrete time, :func:`dlqe` returns the predictor gain :math:`L` by
+default.  For the system :math:`x[n+1] = A x[n] + B u[n] + G w[n]`,
+:math:`y[n] = C x[n] + D u[n] + v[n]`, the prior estimate uses measurements
+through time :math:`n-1` and advances according to
+
+.. math::
+
+   \hat{x}[n+1|n] = A \hat{x}[n|n-1] + B u[n]
+       + L (y[n] - C \hat{x}[n|n-1] - D u[n]).
+
+To update the current state estimate using :math:`y[n]`, set
+`return_filter_form` to True.  This returns the measurement-update gain
+:math:`M = P C^T (C P C^T + RN)^{-1}`, which is related to the predictor
+gain by :math:`L = A M`, and can be computed even when :math:`A` is singular:
+
+.. code::
+
+   M, P, E = ct.dlqe(A, G, C, QN, RN, return_filter_form=True)
+
+The two-step filter uses this gain as follows:
+
+.. math::
+
+   \hat{x}[n|n] &= \hat{x}[n|n-1]
+       + M (y[n] - C \hat{x}[n|n-1] - D u[n]), \\
+   \hat{x}[n+1|n] &= A \hat{x}[n|n] + B u[n].
+
+Both gain forms return the same prior covariance `P` and estimator poles
+`E`.  The posterior covariance is :math:`(I - M C) P`, and the poles are
+the eigenvalues of :math:`A (I - M C)`.  The option is also supported by
+:func:`lqe` when passed a discrete-time system.
 
 .. testsetup:: kalman