Repository: python-control/python-control

Issue:
`solve_flat_ocp` and scalar `timepts`
flowing from gh-1101, which complains about the code in question,

`solve_flat_ocp` argument `timepts` is advertised in the docstring as being allowed to be scalar:

```
    timepts : float or 1D array_like
        The list of points for evaluating cost and constraints, as well as
        the time horizon.  If given as a float, indicates the final time for
        the trajectory (corresponding to xf)
```

but this doesn't work:

```python
import numpy as np
import control as ct
import control.flatsys as fs

f = fs.LinearFlatSystem(ct.ss(ct.tf([1],[1,1])))

def terminal_cost(x, u):
    return (x-5).dot(x-5)+u.dot(u)

traj1 = fs.solve_flat_ocp(f, [0, 1], x0=[23],
                          terminal_cost=terminal_cost)

traj2 = fs.solve_flat_ocp(f, 1, x0=[23],
                          terminal_cost=terminal_cost)

teval = np.linspace(0,1,101)

r1 = traj1.response(teval)
r2 = traj2.response(teval)

assert np.max(abs(r1.x-r2.x)) == 0
assert np.max(abs(r1.u-r2.u)) == 0
assert np.max(abs(r1.y-r2.y)) == 0
```

produces with master @ ebff1259

```
Traceback (most recent call last):
  File "/home/rory/projects/pycontrol/pyflakes/flatbug2.py", line 13, in <module>
    traj2 = fs.solve_flat_ocp(f, 1, x0=[23],
                              terminal_cost=terminal_cost)
  File "/home/rory/src/python-control/control/flatsys/flatsys.py", line 725, in solve_flat_ocp
    T0 = timepts[0] if len(timepts) > 1 else T0
                                             ^^
UnboundLocalError: cannot access local variable 'T0' where it is not associated with a value
```

I think the fix is as below. Diff also removes `Tf`, which is not used in the function.  With this diff the test script runs with error.

```diff
@@ -721,8 +721,7 @@ def solve_flat_ocp(
 
     # Process final time
     timepts = np.atleast_1d(timepts)
-    Tf = timepts[-1]
-    T0 = timepts[0] if len(timepts) > 1 else T0
+    T0 = timepts[0] if len(timepts) > 1 else 0
 
     # Process keyword arguments
     if trajectory_constraints is None:
```

Discussion:

Solution patch:
diff --git a/control/flatsys/flatsys.py b/control/flatsys/flatsys.py
index 7d76b9d78..5818d118b 100644
--- a/control/flatsys/flatsys.py
+++ b/control/flatsys/flatsys.py
@@ -721,8 +721,7 @@ def solve_flat_ocp(
 
     # Process final time
     timepts = np.atleast_1d(timepts)
-    Tf = timepts[-1]
-    T0 = timepts[0] if len(timepts) > 1 else T0
+    T0 = timepts[0] if len(timepts) > 1 else 0
 
     # Process keyword arguments
     if trajectory_constraints is None:
diff --git a/control/xferfcn.py b/control/xferfcn.py
index c9043bd0e..4a8fd4a1c 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1516,6 +1516,12 @@ def _convert_to_transfer_function(
             den = [[[1.] for j in range(sys.ninputs)]
                    for i in range(sys.noutputs)]
         else:
+            # Preallocate numerator and denominator arrays
+            num = [[[] for j in range(sys.ninputs)]
+                   for i in range(sys.noutputs)]
+            den = [[[] for j in range(sys.ninputs)]
+                   for i in range(sys.noutputs)]
+
             try:
                 # Use Slycot to make the transformation
                 # Make sure to convert system matrices to numpy arrays
@@ -1524,12 +1530,6 @@ def _convert_to_transfer_function(
                     sys.nstates, sys.ninputs, sys.noutputs, array(sys.A),
                     array(sys.B), array(sys.C), array(sys.D), tol1=0.0)
 
-                # Preallocate outputs.
-                num = [[[] for j in range(sys.ninputs)]
-                       for i in range(sys.noutputs)]
-                den = [[[] for j in range(sys.ninputs)]
-                       for i in range(sys.noutputs)]
-
                 for i in range(sys.noutputs):
                     for j in range(sys.ninputs):
                         num[i][j] = list(tfout[6][i, j, :])
@@ -1538,16 +1538,13 @@ def _convert_to_transfer_function(
                         den[i][j] = list(tfout[5][i, :])
 
             except ImportError:
-                # If slycot is not available, use signal.lti (SISO only)
-                if sys.ninputs != 1 or sys.noutputs != 1:
-                    raise ControlMIMONotImplemented("Not implemented for " +
-                        "MIMO systems without slycot.")
-
-                # Do the conversion using sp.signal.ss2tf
-                # Note that this returns a 2D array for the numerator
-                num, den = sp.signal.ss2tf(sys.A, sys.B, sys.C, sys.D)
-                num = squeeze(num)  # Convert to 1D array
-                den = squeeze(den)  # Probably not needed
+                # If slycot not available, do conversion using sp.signal.ss2tf
+                for j in range(sys.ninputs):
+                    num_j, den_j = sp.signal.ss2tf(
+                        sys.A, sys.B, sys.C, sys.D, input=j)
+                    for i in range(sys.noutputs):
+                        num[i][j] = num_j[i]
+                        den[i][j] = den_j
 
         newsys = TransferFunction(num, den, sys.dt)
         if use_prefix_suffix:


Test patch:
diff --git a/control/tests/flatsys_test.py b/control/tests/flatsys_test.py
index a12bf1480..5b66edaf5 100644
--- a/control/tests/flatsys_test.py
+++ b/control/tests/flatsys_test.py
@@ -452,6 +452,29 @@ def test_flat_solve_ocp(self, basis):
         np.testing.assert_almost_equal(x_const, x_nlconst)
         np.testing.assert_almost_equal(u_const, u_nlconst)
 
+    def test_solve_flat_ocp_scalar_timepts(self):
+        # scalar timepts gives expected result
+        f = fs.LinearFlatSystem(ct.ss(ct.tf([1],[1,1])))
+
+        def terminal_cost(x, u):
+            return (x-5).dot(x-5)+u.dot(u)
+
+        traj1 = fs.solve_flat_ocp(f, [0, 1], x0=[23],
+                                  terminal_cost=terminal_cost)
+
+        traj2 = fs.solve_flat_ocp(f, 1, x0=[23],
+                                  terminal_cost=terminal_cost)
+
+        teval = np.linspace(0, 1, 101)
+
+        r1 = traj1.response(teval)
+        r2 = traj2.response(teval)
+
+        np.testing.assert_array_equal(r1.x, r2.x)
+        np.testing.assert_array_equal(r1.y, r2.y)
+        np.testing.assert_array_equal(r1.u, r2.u)
+
+
     def test_bezier_basis(self):
         bezier = fs.BezierFamily(4)
         time = np.linspace(0, 1, 100)
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 5c0acd4c8..9cceff389 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -258,7 +258,6 @@ def testNyquist(self, frd_fcn):
         freqplot.nyquist(f1)
         # plt.savefig('/dev/null', format='svg')
 
-    @slycotonly
     @pytest.mark.parametrize(
         "frd_fcn", [ct.frd, ct.FRD, ct.FrequencyResponseData])
     def testMIMO(self, frd_fcn):
@@ -276,7 +275,6 @@ def testMIMO(self, frd_fcn):
             sys.frequency_response(chkpts)[1],
             f1.frequency_response(chkpts)[1])
 
-    @slycotonly
     @pytest.mark.parametrize(
         "frd_fcn", [ct.frd, ct.FRD, ct.FrequencyResponseData])
     def testMIMOfb(self, frd_fcn):
@@ -295,7 +293,6 @@ def testMIMOfb(self, frd_fcn):
             f1.frequency_response(chkpts)[1],
             f2.frequency_response(chkpts)[1])
 
-    @slycotonly
     @pytest.mark.parametrize(
         "frd_fcn", [ct.frd, ct.FRD, ct.FrequencyResponseData])
     def testMIMOfb2(self, frd_fcn):
@@ -316,7 +313,6 @@ def testMIMOfb2(self, frd_fcn):
             f1.frequency_response(chkpts)[1],
             f2.frequency_response(chkpts)[1])
 
-    @slycotonly
     @pytest.mark.parametrize(
         "frd_fcn", [ct.frd, ct.FRD, ct.FrequencyResponseData])
     def testMIMOMult(self, frd_fcn):
@@ -335,7 +331,6 @@ def testMIMOMult(self, frd_fcn):
             (f1*f2).frequency_response(chkpts)[1],
             (sys*sys).frequency_response(chkpts)[1])
 
-    @slycotonly
     @pytest.mark.parametrize(
         "frd_fcn", [ct.frd, ct.FRD, ct.FrequencyResponseData])
     def testMIMOSmooth(self, frd_fcn):
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 78177995d..48aa71795 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2288,7 +2288,6 @@ def test_signal_indexing():
         resp.outputs['y[0]', 'u[0]']
 
 
-@slycotonly
 @pytest.mark.parametrize("fcn, spec, expected, missing", [
     (ct.ss, {}, "states=4, outputs=3, inputs=2", r"dt|name"),
     (ct.tf, {}, "outputs=3, inputs=2", r"dt|states|name"),
diff --git a/control/tests/lti_test.py b/control/tests/lti_test.py
index e93138af3..1b059ef6f 100644
--- a/control/tests/lti_test.py
+++ b/control/tests/lti_test.py
@@ -309,7 +309,6 @@ def test_squeeze_exceptions(self, fcn):
             evalfr(sys, [[0.1j, 1j], [1j, 10j]])
 
 
-@slycotonly
 @pytest.mark.parametrize(
     "outdx, inpdx, key",
     [('y[0]', 'u[1]', (0, 1)),
@@ -356,7 +355,6 @@ def test_subsys_indexing(fcn, outdx, inpdx, key):
                 subsys_chk.frequency_response(omega).response)
 
 
-@slycotonly
 @pytest.mark.parametrize("op", [
     '__mul__', '__rmul__', '__add__', '__radd__', '__sub__', '__rsub__'])
 @pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.frd])
diff --git a/control/tests/matlab2_test.py b/control/tests/matlab2_test.py
index 5eedfc2ec..d7c0c1ce1 100644
--- a/control/tests/matlab2_test.py
+++ b/control/tests/matlab2_test.py
@@ -49,7 +49,6 @@ def MIMO_mats(self):
         D = zeros((2, 2))
         return A, B, C, D
 
-    @slycotonly
     def test_dcgain_mimo(self, MIMO_mats):
         """Test function dcgain with MIMO systems"""
         #Test MIMO systems
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 366c5ef3b..6d57a38a9 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -248,7 +248,6 @@ def test_zero_siso(self, sys222):
 
         np.testing.assert_almost_equal(true_z, z)
 
-    @slycotonly
     def test_zero_mimo_sys322_square(self, sys322):
         """Evaluate the zeros of a square MIMO system."""
 
@@ -256,7 +255,6 @@ def test_zero_mimo_sys322_square(self, sys322):
         true_z = np.sort([44.41465, -0.490252, -5.924398])
         np.testing.assert_array_almost_equal(z, true_z)
 
-    @slycotonly
     def test_zero_mimo_sys222_square(self, sys222):
         """Evaluate the zeros of a square MIMO system."""
 
@@ -320,7 +318,6 @@ def test_multiply_ss(self, sys222, sys322):
         np.testing.assert_array_almost_equal(sys.C, C)
         np.testing.assert_array_almost_equal(sys.D, D)
 
-    @slycotonly
     def test_add_sub_mimo_siso(self):
         # Test SS with SS
         ss_siso = StateSpace(
@@ -702,8 +699,6 @@ def test_call(self, dt, omega, resp):
         with pytest.raises(AttributeError):
             sys.evalfr(omega)
 
-
-    @slycotonly
     def test_freq_resp(self):
         """Evaluate the frequency response at multiple frequencies."""
 
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index cab556cc2..00000f52c 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -10,9 +10,9 @@
 import pytest
 
 import control as ct
-from control import (StateSpace, TransferFunction, defaults, evalfr, isctime,
-                     isdtime, reset_defaults, rss, sample_system, set_defaults,
-                     ss, ss2tf, tf, tf2ss, zpk)
+from control import StateSpace, TransferFunction, defaults, evalfr, isctime, \
+    isdtime, reset_defaults, rss, sample_system, set_defaults, ss, ss2tf, tf, \
+    tf2ss, zpk
 from control.statesp import _convert_to_statespace
 from control.tests.conftest import slycotonly
 from control.xferfcn import _convert_to_transfer_function, _tf_close_coeff
@@ -186,7 +186,6 @@ def test_reverse_sign_siso(self):
         np.testing.assert_allclose(sys2.num, [[[-1., -3., -5.]]])
         np.testing.assert_allclose(sys2.den, [[[1., 6., 2., -1.]]])
 
-    @slycotonly
     def test_reverse_sign_mimo(self):
         """Negate a MIMO system."""
         num1 = [[[1., 2.], [0., 3.], [2., -1.]],
@@ -228,7 +227,6 @@ def test_add_siso(self):
         np.testing.assert_allclose(sys3.num, [[[20., 4., -8]]])
         np.testing.assert_allclose(sys3.den, [[[1., 6., 1., -7., -2., 1.]]])
 
-    @slycotonly
     def test_add_mimo(self):
         """Add two MIMO systems."""
         num1 = [[[1., 2.], [0., 3.], [2., -1.]],
@@ -276,7 +274,6 @@ def test_subtract_siso(self):
         np.testing.assert_allclose(sys4.num, [[[-2., -6., 12., 10., 2.]]])
         np.testing.assert_allclose(sys4.den, [[[1., 6., 1., -7., -2., 1.]]])
 
-    @slycotonly
     def test_subtract_mimo(self):
         """Subtract two MIMO systems."""
         num1 = [[[1., 2.], [0., 3.], [2., -1.]],
@@ -327,7 +324,6 @@ def test_multiply_siso(self):
         np.testing.assert_allclose(sys3.num, sys4.num)
         np.testing.assert_allclose(sys3.den, sys4.den)
 
-    @slycotonly
     def test_multiply_mimo(self):
         """Multiply two MIMO systems."""
         num1 = [[[1., 2.], [0., 3.], [2., -1.]],
@@ -714,7 +710,6 @@ def test_call_dtime(self):
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1], 0.1)
         np.testing.assert_array_almost_equal(sys(1j), -0.5 - 0.5j)
 
-    @slycotonly
     def test_call_mimo(self):
         """Evaluate the frequency response of a MIMO system at one frequency."""
 
@@ -755,7 +750,6 @@ def test_frequency_response_siso(self):
         np.testing.assert_array_almost_equal(phase, truephase)
         np.testing.assert_array_almost_equal(omega, trueomega)
 
-    @slycotonly
     def test_freqresp_mimo(self):
         """Evaluate the MIMO magnitude and phase at multiple frequencies."""
         num = [[[1., 2.], [0., 3.], [2., -1.]],
@@ -852,7 +846,6 @@ def test_common_den_nonproper(self):
         _, den2, _ = tf2._common_den(allow_nonproper=True)
         np.testing.assert_array_almost_equal(den2, common_den_ref)
 
-    @slycotonly
     def test_pole_mimo(self):
         """Test for correct MIMO poles."""
         sys = TransferFunction(
@@ -936,7 +929,6 @@ def test_append(self):
         tf_appended_2 = tf1.append(tf2).append(tf3)
         assert _tf_close_coeff(tf_exp_2, tf_appended_2)
 
-    @slycotonly
     def test_convert_to_transfer_function(self):
         """Test for correct state space to transfer function conversion."""
         A = [[1., -2.], [-3., 4.]]
@@ -1023,7 +1015,6 @@ def test_state_space_conversion_mimo(self):
         np.testing.assert_array_almost_equal(H.num[1][0], H2.num[1][0])
         np.testing.assert_array_almost_equal(H.den[1][0], H2.den[1][0])
 
-    @slycotonly
     def test_indexing(self):
         """Test TF scalar indexing and slice"""
         tm = ss2tf(rss(5, 3, 3))
@@ -1213,7 +1204,6 @@ def test_printing_polynomial(self, args, outputfmt, var, dt, dtstring):
             assert len(polystr[0].split('\n')) == 4
         assert polystr[2] == outputfmt.format(var=var)
 
-    @slycotonly
     def test_printing_mimo(self):
         """Print MIMO, continuous time"""
         sys = ss2tf(rss(4, 2, 3))
@@ -1332,7 +1322,6 @@ def test_printing_zpk_mimo(self, num, den, output):
         res = str(G)
         assert res.partition('\n\n')[2] == output
 
-    @slycotonly
     def test_size_mismatch(self):
         """Test size mismacht"""
         sys1 = ss2tf(rss(2, 2, 2))