Repository: python-control/python-control
Pull request #1125 (2025-02-09 05:48:33)

Uniform processing of time response and optimization parameters
This PR regularizes the use of various keywords for time responses and optimization routines, while maintaining backward compatibility with existing code.  The motivation for the PR is that over time we have introduced different terms to refer to parameters that have the same purpose.  This can be confusing for people trying to use the package.

For example:
* The named parameters `x0` and `X0` are used in different places for the initial state in simulation and optimization functions.   
* The parameter `return_x` is used when a function should return the state as part of a tuple, though sometime `return_states` is used instead.  
* Responses use the `states` property (not `x`) and setting the the state names is done via `states`. 
* Similar issues arise for inputs (`u0`, `U`,  `inputs`) and outputs (`y0`, `outputs`).
* Time points for simulation and optimization are sometimes referred to as `timepts` but other times as `T`.
* In the optimization and flatness modules, cost functions are referred to as `cost`, `trajectory_cost` or `integral cost`, depending on the function, and constraints are referred to as `trajectory_constraints` in some places as `constraints` in others.
* Other parameters have names that don't seem to fit existing python-control naming patterns: `t_eval`, `yfinal`, etc.

To address these issues, this PR introduces a common set of terms across various time response and optimal control functions, with the following more consistent usage patterns and functionality:

* Parameters specifying the inputs, outputs, and states are referred to as `inputs`, `outputs`, and `states` consistently throughout the functions.
* Variables associated with inputs, outputs, states and time use those words plus an appropriate modifier: `initial_state`, `final_output`, `input_indices` etc.
* Aliases are used both to maintain backward compatibility and to allow shorthand descriptions: e.g., `U`, `Y`, `X0`.  Short form aliases are documented in docstrings by listing the parameter as `long_form (or sf) : type`.
* Existing legacy keywords are allowed and generate a `PendingDeprecationWarning`.
* Specifying a parameter value in two different ways (eg, via long form and an alias) generates a `TypeError`.

The implementation is done in a manner that can later be utilized for other functions where we want to have aliases and legacy keys, as documented in the Developer Notes section of the Reference Manual:

> [P]arameter names are generally longer strings that describe the purpose fo the paramater. Similar to `matplotlib` (e.g., the use of `lw` as an alias for `linewidth`), some commonly used parameter names can be specified using an “alias” that allows the use of a shorter key.
>
>Named parameter and keyword variable aliases are processed using the `config._process_kwargs()`  and `config._process_param()` functions. These functions allow the specification of a list of aliases and a list of legacy keys for a given named parameter or keyword. To make use of these functions, the `_process_kwargs()` is first called to update the kwargs variable by replacing aliases with the full key:
> ```
> _process_kwargs(kwargs, aliases)
> ```
>The values for named parameters can then be assigned to a local variable using a call to process_param() of the form:
> ```
> var = _process_kwargs('param', param, kwargs, aliases)
> ```
> where `param` is the named parameter used in the function signature and `var` is the local variable in the function (may also be param, but doesn’t have to be).
>
> ...
> 
> The alias mapping is a dictionary that returns a tuple consisting of valid aliases and legacy aliases:
> ```
> alias_mapping = {
>     'argument_name_1', (['alias', ...], ['legacy', ...]),
>      ...}
> ```
> If an alias is present in the dictionary of keywords, it will be used to set the value of the argument. If a legacy keyword is used, a warning is issued.

Two primary tables have been created in this PR.  The first is a 
dictionary of aliases and legacy names for time response commands (from `timeresp.py`):
```
_timeresp_aliases = {
    # param:            ([alias, ...], [legacy, ...])
    'timepts':          (['T'],        []),
    'inputs':           (['U'],        ['u']),
    'outputs':          (['Y'],        ['y']),
    'initial_state':    (['X0'],       ['x0']),
    'final_output':     (['yfinal'],   []),
    'return_states':    (['return_x'], []),
    'evaluation_times': (['t_eval'],   []),
    'timepts_num':      (['T_num'],    []),
    'input_indices':    (['input'],    []),
    'output_indices':   (['output'],   []),
}
```
The second is a dictionary of aliases for optimal control functions (from `optimal.py`):
```
_optimal_aliases = {
    # param:                  ([alias, ...],                [legacy, ...])
    'integral_cost':          (['trajectory_cost', 'cost'], []),
    'initial_state':          (['x0', 'X0'],                []),
    'initial_input':          (['u0', 'U0'],                []),
    'final_state':            (['xf'],                      []),
    'final_input':            (['uf'],                      []),
    'initial_time':           (['T0'],                      []),
    'trajectory_constraints': (['constraints'],             []),
    'return_states':          (['return_x'],                []),
}
```



Patch:
diff --git a/control/config.py b/control/config.py
index 85284912..8da7e2fc 100644
--- a/control/config.py
+++ b/control/config.py
@@ -1,7 +1,7 @@
 # config.py - package defaults
 # RMM, 4 Nov 2012
 #
-# TODO: add ability to read/write configuration files (ala matplotlib)
+# TODO: add ability to read/write configuration files (a la matplotlib)
 
 """Functions to access default parameter values.
 
@@ -384,10 +384,13 @@ def use_legacy_defaults(version):
 def _process_legacy_keyword(kwargs, oldkey, newkey, newval, warn_oldkey=True):
     """Utility function for processing legacy keywords.
 
+    .. deprecated:: 0.10.2
+        Replace with `_process_param` or `_process_kwargs`.
+
     Use this function to handle a legacy keyword that has been renamed.
     This function pops the old keyword off of the kwargs dictionary and
     issues a warning.  If both the old and new keyword are present, a
-    ControlArgument exception is raised.
+    `ControlArgument` exception is raised.
 
     Parameters
     ----------
@@ -412,6 +415,10 @@ def _process_legacy_keyword(kwargs, oldkey, newkey, newval, warn_oldkey=True):
         Value of the (new) keyword.
 
     """
+    # TODO: turn on this warning when ready to deprecate
+    # warnings.warn(
+    #     "replace `_process_legacy_keyword` with `_process_param` "
+    #     "or `_process_kwargs`", PendingDeprecationWarning)
     if oldkey in kwargs:
         if warn_oldkey:
             warnings.warn(
@@ -424,3 +431,143 @@ def _process_legacy_keyword(kwargs, oldkey, newkey, newval, warn_oldkey=True):
             return kwargs.pop(oldkey)
     else:
         return newval
+
+
+def _process_param(name, defval, kwargs, alias_mapping, sigval=None):
+    """Process named parameter, checking aliases and legacy usage.
+
+    Helper function to process function arguments by mapping aliases to
+    either their default keywords or to a named argument.  The alias
+    mapping is a dictionary that returns a tuple consisting of valid
+    aliases and legacy aliases::
+
+       alias_mapping = {
+            'argument_name_1': (['alias', ...], ['legacy', ...]),
+            ...}
+
+    If `param` is a named keyword in the function signature with default
+    value `defval`, a typical calling sequence at the start of a function
+    is::
+
+        param = _process_param('param', defval, kwargs, function_aliases)
+
+    If `param` is a variable keyword argument (in `kwargs`), `defval` can
+    be passed as either None or the default value to use if `param` is not
+    present in `kwargs`.
+
+    Parameters
+    ----------
+    name : str
+        Name of the parameter to be checked.
+    defval : object or dict
+        Default value for the parameter.
+    kwargs : dict
+        Dictionary of variable keyword arguments.
+    alias_mapping : dict
+        Dictionary providing aliases and legacy names.
+    sigval : object, optional
+        Default value specified in the function signature (default = None).
+        If specified, an error will be generated if `defval` is different
+        than `sigval` and an alias or legacy keyword is given.
+
+    Returns
+    -------
+    newval : object
+        New value of the named parameter.
+
+    Raises
+    ------
+    TypeError
+        If multiple keyword aliases are used for the same parameter.
+
+    Warns
+    -----
+    PendingDeprecationWarning
+        If legacy name is used to set the value for the variable.
+
+    """
+    # Check to see if the parameter is in the keyword list
+    if name in kwargs:
+        if defval != sigval:
+            raise TypeError(f"multiple values for parameter {name}")
+        newval = kwargs.pop(name)
+    else:
+        newval = defval
+
+    # Get the list of aliases and legacy names
+    aliases, legacy = alias_mapping[name]
+
+    for kw in legacy:
+        if kw in kwargs:
+            warnings.warn(
+                f"alias `{kw}` is legacy name; use `{name}` instead",
+                PendingDeprecationWarning)
+            kwval = kwargs.pop(kw)
+            if newval != defval and kwval != newval:
+                raise TypeError(
+                    f"multiple values for parameter `{name}` (via {kw})")
+            newval = kwval
+
+    for kw in aliases:
+        if kw in kwargs:
+            kwval = kwargs.pop(kw)
+            if newval != defval and kwval != newval:
+                raise TypeError(
+                    f"multiple values for parameter `{name}` (via {kw})")
+            newval = kwval
+
+    return newval
+
+
+def _process_kwargs(kwargs, alias_mapping):
+    """Process aliases and legacy keywords.
+
+    Helper function to process function arguments by mapping aliases to
+    their default keywords.  The alias mapping is a dictionary that returns
+    a tuple consisting of valid aliases and legacy aliases::
+
+       alias_mapping = {
+            'argument_name_1': (['alias', ...], ['legacy', ...]),
+            ...}
+
+    If an alias is present in the dictionary of keywords, it will be used
+    to set the value of the argument.  If a legacy keyword is used, a
+    warning is issued.
+
+    Parameters
+    ----------
+    kwargs : dict
+        Dictionary of variable keyword arguments.
+    alias_mapping : dict
+        Dictionary providing aliases and legacy names.
+
+    Raises
+    ------
+    TypeError
+        If multiple keyword aliased are used for the same parameter.
+
+    Warns
+    -----
+    PendingDeprecationWarning
+        If legacy name is used to set the value for the variable.
+
+    """
+    for name in alias_mapping or []:
+        aliases, legacy = alias_mapping[name]
+
+        for kw in legacy:
+            if kw in kwargs:
+                warnings.warn(
+                    f"alias `{kw}` is legacy name; use `{name}` instead",
+                    PendingDeprecationWarning)
+                if name in kwargs:
+                    raise TypeError(
+                        f"multiple values for parameter `{name}` (via {kw})")
+                kwargs[name] = kwargs.pop(kw)
+
+        for kw in aliases:
+            if kw in kwargs:
+                if name in kwargs:
+                    raise TypeError(
+                        f"multiple values for parameter `{name}` (via {kw})")
+                kwargs[name] = kwargs.pop(kw)
diff --git a/control/flatsys/flatsys.py b/control/flatsys/flatsys.py
index b57f9bd7..92d32d01 100644
--- a/control/flatsys/flatsys.py
+++ b/control/flatsys/flatsys.py
@@ -12,9 +12,10 @@ import numpy as np
 import scipy as sp
 import scipy.optimize
 
-from ..config import _process_legacy_keyword
+from ..config import _process_kwargs, _process_param
 from ..exception import ControlArgument
 from ..nlsys import NonlinearIOSystem
+from ..optimal import _optimal_aliases
 from ..timeresp import _check_convert_array
 from .poly import PolyFamily
 from .systraj import SystemTrajectory
@@ -325,7 +326,8 @@ l
 
 # Solve a point to point trajectory generation problem for a flat system
 def point_to_point(
-        sys, timepts, x0=0, u0=0, xf=0, uf=0, T0=0, cost=None, basis=None,
+        sys, timepts, initial_state=0, initial_input=0, final_state=0,
+        final_input=0, initial_time=0, integral_cost=None, basis=None,
         trajectory_constraints=None, initial_guess=None, params=None, **kwargs):
     """Compute trajectory between an initial and final conditions.
 
@@ -340,32 +342,30 @@ def point_to_point(
         and produces the flag of flat outputs and a function
         `~FlatSystem.reverse` that takes the flag of the flat output and
         produces the state and input.
-
     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)
-
-    x0, u0, xf, uf : 1D arrays
-        Define the desired initial and final conditions for the system.  If
-        any of the values are given as None, they are replaced by a vector of
-        zeros of the appropriate dimension.
-
-    T0 : float, optional
+    initial_state (or x0) : 1D array_like
+        Initial state for the system.  Defaults to zero.
+    initial_input (or u0) : 1D array_like
+        Initial input for the system.  Defaults to zero.
+    final_state (or xf) : 1D array_like
+        Final state for the system.  Defaults to zero.
+    final_input (or uf) : 1D array_like
+        Final input for the system.  Defaults to zero.
+    initial_time (or T0) : float, optional
         The initial time for the trajectory (corresponding to x0).  If not
         specified, its value is taken to be zero.
-
     basis : `BasisFamily` object, optional
         The basis functions to use for generating the trajectory.  If not
         specified, the `PolyFamily` basis family
         will be used, with the minimal number of elements required to find a
         feasible trajectory (twice the number of system states)
-
-    cost : callable
+    integral_cost (or cost) : callable
         Function that returns the integral cost given the current state
-        and input.  Called as ``cost(x, u)``.
-
-    trajectory_constraints : list of tuples, optional
+        and input.  Called as ``integral_cost(x, u)``.
+    trajectory_constraints (or constraints) : list of tuples, optional
         List of constraints that should hold at each point in the time
         vector.  Each element of the list should consist of a tuple with
         first element given by `scipy.optimize.LinearConstraint` or
@@ -382,24 +382,13 @@ def point_to_point(
           trajectory and compared against the upper and lower bounds.
 
         The constraints are applied at each time point along the trajectory.
-
     initial_guess : 2D array_like, optional
         Initial guess for the trajectory coefficients (not implemented).
-
     params : dict, optional
         Parameter values for the system.  Passed to the evaluation
         functions for the system as default values, overriding internal
         defaults.
 
-    minimize_method : str, optional
-        Set the method used by `scipy.optimize.minimize`.
-
-    minimize_options : str, optional
-        Set the options keyword used by `scipy.optimize.minimize`.
-
-    minimize_kwargs : str, optional
-        Pass additional keywords to `scipy.optimize.minimize`.
-
     Returns
     -------
     traj : `SystemTrajectory` object
@@ -407,6 +396,15 @@ def point_to_point(
         `~SystemTrajectory.eval` function, we can be used to
         compute the value of the state and input and a given time t.
 
+    Other Parameters
+    ----------------
+    minimize_method : str, optional
+        Set the method used by `scipy.optimize.minimize`.
+    minimize_options : str, optional
+        Set the options keyword used by `scipy.optimize.minimize`.
+    minimize_kwargs : str, optional
+        Pass additional keywords to `scipy.optimize.minimize`.
+
     Notes
     -----
     Additional keyword parameters can be used to fine tune the behavior of
@@ -414,6 +412,24 @@ def point_to_point(
     `OptimalControlProblem` for more information.
 
     """
+    # Process parameter and keyword arguments
+    _process_kwargs(kwargs, _optimal_aliases)
+    x0 = _process_param(
+        'initial_state', initial_state, kwargs, _optimal_aliases, sigval=0)
+    u0 = _process_param(
+        'initial_input', initial_input, kwargs, _optimal_aliases, sigval=0)
+    xf = _process_param(
+        'final_state', final_state, kwargs, _optimal_aliases, sigval=0)
+    uf = _process_param(
+        'final_input', final_input, kwargs, _optimal_aliases, sigval=0)
+    T0 = _process_param(
+        'initial_time', initial_time, kwargs, _optimal_aliases, sigval=0)
+    cost = _process_param(
+        'integral_cost', integral_cost, kwargs, _optimal_aliases)
+    trajectory_constraints = _process_param(
+        'trajectory_constraints', trajectory_constraints, kwargs,
+        _optimal_aliases)
+
     #
     # Make sure the problem is one that we can handle
     #
@@ -431,13 +447,6 @@ def point_to_point(
     Tf = timepts[-1]
     T0 = timepts[0] if len(timepts) > 1 else T0
 
-    # Process keyword arguments
-    trajectory_constraints = _process_legacy_keyword(
-        kwargs, 'constraints', 'trajectory_constraints',
-        trajectory_constraints, warn_oldkey=False)
-    cost = _process_legacy_keyword(
-        kwargs, 'trajectory_cost', 'cost', cost, warn_oldkey=False)
-
     minimize_kwargs = {}
     minimize_kwargs['method'] = kwargs.pop('minimize_method', None)
     minimize_kwargs['options'] = kwargs.pop('minimize_options', {})
@@ -657,8 +666,8 @@ def point_to_point(
 
 # Solve a point to point trajectory generation problem for a flat system
 def solve_flat_optimal(
-        sys, timepts, x0=0, u0=0, trajectory_cost=None, basis=None,
-        terminal_cost=None, trajectory_constraints=None,
+        sys, timepts, initial_state=0, initial_input=0, integral_cost=None,
+        basis=None, terminal_cost=None, trajectory_constraints=None,
         initial_guess=None, params=None, **kwargs):
     """Compute trajectory between an initial and final conditions.
 
@@ -673,31 +682,25 @@ def solve_flat_optimal(
         and produces the flag of flat outputs and a function
         `~FlatSystem.reverse` that takes the flag of the flat output and
         produces the state and input.
-
     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)
-
-    x0, u0 : 1D arrays
-        Define the initial conditions for the system.  If either of the
-        values are given as None, they are replaced by a vector of zeros of
-        the appropriate dimension.
-
+    initial_state (or x0), input_input (or u0) : 1D arrays
+        Define the initial conditions for the system (default = 0).
+    initial_input (or u0) : 1D array_like
+        Initial input for the system.  Defaults to zero.
     basis : `BasisFamily` object, optional
         The basis functions to use for generating the trajectory.  If not
         specified, the `PolyFamily` basis family
         will be used, with the minimal number of elements required to find a
         feasible trajectory (twice the number of system states)
-
-    trajectory_cost : callable
+    integral_cost : callable
         Function that returns the integral cost given the current state
         and input.  Called as ``cost(x, u)``.
-
     terminal_cost : callable
         Function that returns the terminal cost given the state and input.
         Called as ``cost(x, u)``.
-
     trajectory_constraints : list of tuples, optional
         List of constraints that should hold at each point in the time
         vector.  Each element of the list should consist of a tuple with
@@ -715,15 +718,22 @@ def solve_flat_optimal(
           trajectory and compared against the upper and lower bounds.
 
         The constraints are applied at each time point along the trajectory.
-
     initial_guess : 2D array_like, optional
         Initial guess for the optimal trajectory of the flat outputs.
-
     params : dict, optional
         Parameter values for the system.  Passed to the evaluation
         functions for the system as default values, overriding internal
         defaults.
 
+    Returns
+    -------
+    traj : `SystemTrajectory`
+        The system trajectory is returned as an object that implements the
+        `SystemTrajectory.eval` function, we can be used to
+        compute the value of the state and input and a given time `t`.
+
+    Other Parameters
+    ----------------
     minimize_method : str, optional
         Set the method used by `scipy.optimize.minimize`.
 
@@ -733,13 +743,6 @@ def solve_flat_optimal(
     minimize_kwargs : str, optional
         Pass additional keywords to `scipy.optimize.minimize`.
 
-    Returns
-    -------
-    traj : `SystemTrajectory`
-        The system trajectory is returned as an object that implements the
-        `SystemTrajectory.eval` function, we can be used to
-        compute the value of the state and input and a given time `t`.
-
     Notes
     -----
     Additional keyword parameters can be used to fine tune the behavior of
@@ -758,6 +761,18 @@ def solve_flat_optimal(
     used to overcome these errors.
 
     """
+    # Process parameter and keyword arguments
+    _process_kwargs(kwargs, _optimal_aliases)
+    x0 = _process_param(
+        'initial_state', initial_state, kwargs, _optimal_aliases, sigval=0)
+    u0 = _process_param(
+        'initial_input', initial_input, kwargs, _optimal_aliases, sigval=0)
+    trajectory_cost = _process_param(
+        'integral_cost', integral_cost, kwargs, _optimal_aliases)
+    trajectory_constraints = _process_param(
+        'trajectory_constraints', trajectory_constraints, kwargs,
+        _optimal_aliases)
+
     #
     # Make sure the problem is one that we can handle
     #
@@ -770,12 +785,6 @@ def solve_flat_optimal(
     timepts = np.atleast_1d(timepts)
     T0 = timepts[0] if len(timepts) > 1 else 0
 
-    # Process keyword arguments
-    trajectory_constraints = _process_legacy_keyword(
-        kwargs, 'constraints', 'trajectory_constraints', trajectory_constraints)
-    trajectory_cost = _process_legacy_keyword(
-        kwargs, 'cost', 'trajectory_cost', trajectory_cost)
-
     minimize_kwargs = {}
     minimize_kwargs['method'] = kwargs.pop('minimize_method', None)
     minimize_kwargs['options'] = kwargs.pop('minimize_options', {})
diff --git a/control/nlsys.py b/control/nlsys.py
index 32524a9c..30f06f81 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -22,10 +22,11 @@ import numpy as np
 import scipy as sp
 
 from . import config
+from .config import _process_param, _process_kwargs
 from .iosys import InputOutputSystem, _parse_spec, _process_iosys_keywords, \
     common_timebase, iosys_repr, isctime, isdtime
 from .timeresp import TimeResponseData, TimeResponseList, \
-    _check_convert_array, _process_time_response
+    _check_convert_array, _process_time_response, _timeresp_aliases
 
 __all__ = ['NonlinearIOSystem', 'InterconnectedSystem', 'nlsys',
            'input_output_response', 'find_eqpt', 'linearize',
@@ -1471,9 +1472,9 @@ def nlsys(updfcn, outfcn=None, **kwargs):
 
 
 def input_output_response(
-        sys, T, U=0., X0=0, params=None, ignore_errors=False,
-        transpose=False, return_x=False, squeeze=None,
-        solve_ivp_kwargs=None, t_eval='T', **kwargs):
+        sys, timepts=None, inputs=0., initial_state=0., params=None,
+        ignore_errors=False, transpose=False, return_states=False,
+        squeeze=None, solve_ivp_kwargs=None, evaluation_times='T', **kwargs):
     """Compute the output response of a system to a given input.
 
     Simulate a dynamical system with a given input and return its output
@@ -1483,22 +1484,23 @@ def input_output_response(
     ----------
     sys : `NonlinearIOSystem` or list of `NonlinearIOSystem`
         I/O system(s) for which input/output response is simulated.
-    T : array_like
+    timepts (or T) : array_like
         Time steps at which the input is defined; values must be evenly spaced.
-    U : array_like, list, or number, optional
-        Input array giving input at each time `T` (default = 0).  If a list
-        is specified, each element in the list will be treated as a portion
-        of the input and broadcast (if necessary) to match the time vector.
-    X0 : array_like, list, or number, optional
+    inputs (or U) : array_like, list, or number, optional
+        Input array giving input at each time in `timepts` (default =
+        0). If a list is specified, each element in the list will be
+        treated as a portion of the input and broadcast (if necessary) to
+        match the time vector.
+    initial_state (or X0) : array_like, list, or number, optional
         Initial condition (default = 0).  If a list is given, each element
         in the list will be flattened and stacked into the initial
         condition.  If a smaller number of elements are given that the
         number of states in the system, the initial condition will be padded
         with zeros.
-    t_eval : array-list, optional
+    evaluation_times (or t_eval) : array-list, optional
         List of times at which the time response should be computed.
-        Defaults to `T`.
-    return_x : bool, optional
+        Defaults to `timepts`.
+    return_states (or return_x) : bool, optional
         If True, return the state vector when assigning to a tuple.  See
         `forced_response` for more details.  If True, return the values of
         the state at each time Default is False.
@@ -1523,7 +1525,7 @@ def input_output_response(
         method.  See `TimeResponseData` for more detailed information.
     response.time : array
         Time values of the output.
-    response.output : array
+    response.outputs : array
         Response of the system.  If the system is SISO and `squeeze` is not
         True, the array is 1D (indexed by time).  If the system is not SISO
         or `squeeze` is False, the array is 2D (indexed by output and time).
@@ -1581,6 +1583,18 @@ def input_output_response(
     #
     # Process keyword arguments
     #
+    _process_kwargs(kwargs, _timeresp_aliases)
+    T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+    U = _process_param('inputs', inputs, kwargs, _timeresp_aliases, sigval=0.)
+    X0 = _process_param(
+        'initial_state', initial_state, kwargs, _timeresp_aliases, sigval=0.)
+    return_x = _process_param(
+        'return_states', return_states, kwargs, _timeresp_aliases,
+        sigval=False)
+    # TODO: replace default value of evaluation_times with None?
+    t_eval = _process_param(
+        'evaluation_times', evaluation_times, kwargs, _timeresp_aliases,
+        sigval='T')
 
     # Figure out the method to be used
     solve_ivp_kwargs = solve_ivp_kwargs.copy() if solve_ivp_kwargs else {}
@@ -1605,9 +1619,10 @@ def input_output_response(
         sysdata, responses = sys, []
         for sys in sysdata:
             responses.append(input_output_response(
-                sys, T, U=U, X0=X0, params=params, transpose=transpose,
-                return_x=return_x, squeeze=squeeze, t_eval=t_eval,
-                solve_ivp_kwargs=solve_ivp_kwargs, **kwargs))
+                sys, timepts=T, inputs=U, initial_state=X0, params=params,
+                transpose=transpose, return_states=return_x, squeeze=squeeze,
+                evaluation_times=t_eval, solve_ivp_kwargs=solve_ivp_kwargs,
+                **kwargs))
         return TimeResponseList(responses)
 
     # Sanity checking on the input
@@ -1894,8 +1909,9 @@ class OperatingPoint():
 
 
 def find_operating_point(
-        sys, x0, u0=None, y0=None, t=0, params=None, iu=None, iy=None,
-        ix=None, idx=None, dx0=None, root_method=None, root_kwargs=None,
+        sys, initial_state=0., inputs=None, outputs=None, t=0, params=None,
+        input_indices=None, output_indices=None, state_indices=None,
+        deriv_indices=None, derivs=None, root_method=None, root_kwargs=None,
         return_outputs=None, return_result=None, **kwargs):
     """Find an operating point for an input/output system.
 
@@ -1929,13 +1945,13 @@ def find_operating_point(
     ----------
     sys : `NonlinearIOSystem`
         I/O system for which the operating point is sought.
-    x0 : list of initial state values
+    initial_state (or x0) : list of initial state values
         Initial guess for the value of the state near the operating point.
-    u0 : list of input values, optional
+    inputs (or u0) : list of input values, optional
         If `y0` is not specified, sets the value of the input.  If `y0` is
         given, provides an initial guess for the value of the input.  Can
         be omitted if the system does not have any inputs.
-    y0 : list of output values, optional
+    outputs (or y0) : list of output values, optional
         If specified, sets the desired values of the outputs at the
         operating point.
     t : float, optional
@@ -1943,22 +1959,22 @@ def find_operating_point(
     params : dict, optional
         Parameter values for the system.  Passed to the evaluation functions
         for the system as default values, overriding internal defaults.
-    iu : list of input indices, optional
+    input_indices (or iu) : list of input indices, optional
         If specified, only the inputs with the given indices will be fixed at
         the specified values in solving for an operating point.  All other
         inputs will be varied.  Input indices can be listed in any order.
-    iy : list of output indices, optional
+    output_indices (or iy) : list of output indices, optional
         If specified, only the outputs with the given indices will be fixed
         at the specified values in solving for an operating point.  All other
         outputs will be varied.  Output indices can be listed in any order.
-    ix : list of state indices, optional
+    state_indices (or ix) : list of state indices, optional
         If specified, states with the given indices will be fixed at the
         specified values in solving for an operating point.  All other
         states will be varied.  State indices can be listed in any order.
-    dx0 : list of update values, optional
+    derivs (or dx0) : list of update values, optional
         If specified, the value of update map must match the listed value
         instead of the default value for an equilibrium point.
-    idx : list of state indices, optional
+    deriv_indices (or idx) : list of state indices, optional
         If specified, state updates with the given indices will have their
         update maps fixed at the values given in `dx0`.  All other update
         values will be ignored in solving for an operating point.  State
@@ -2014,8 +2030,29 @@ def find_operating_point(
     from scipy.optimize import root
 
     # Process keyword arguments
-    return_outputs = config._process_legacy_keyword(
-        kwargs, 'return_y', 'return_outputs', return_outputs)
+    aliases = {
+        'initial_state': (['x0', 'X0'], []),
+        'inputs': (['u0'], []),
+        'outputs': (['y0'], []),
+        'derivs': (['dx0'], []),
+        'input_indices': (['iu'], []),
+        'output_indices': (['iy'], []),
+        'state_indices': (['ix'], []),
+        'deriv_indices': (['idx'], []),
+        'return_outputs': ([], ['return_y']),
+    }
+    _process_kwargs(kwargs, aliases)
+    x0 = _process_param(
+        'initial_state', initial_state, kwargs, aliases, sigval=0.)
+    u0 = _process_param('inputs', inputs, kwargs, aliases)
+    y0 = _process_param('outputs', outputs, kwargs, aliases)
+    dx0 = _process_param('derivs', derivs, kwargs, aliases)
+    iu = _process_param('input_indices', input_indices, kwargs, aliases)
+    iy = _process_param('output_indices', output_indices, kwargs, aliases)
+    ix = _process_param('state_indices', state_indices, kwargs, aliases)
+    idx = _process_param('deriv_indices', deriv_indices, kwargs, aliases)
+    return_outputs = _process_param(
+        'return_outputs', return_outputs, kwargs, aliases)
     if kwargs:
         raise TypeError("unrecognized keyword(s): " + str(kwargs))
 
@@ -2025,9 +2062,9 @@ def find_operating_point(
         root_kwargs['method'] = root_method
 
     # Figure out the number of states, inputs, and outputs
-    x0, nstates = _process_vector_argument(x0, "x0", sys.nstates)
-    u0, ninputs = _process_vector_argument(u0, "u0", sys.ninputs)
-    y0, noutputs = _process_vector_argument(y0, "y0", sys.noutputs)
+    x0, nstates = _process_vector_argument(x0, "initial_states", sys.nstates)
+    u0, ninputs = _process_vector_argument(u0, "inputs", sys.ninputs)
+    y0, noutputs = _process_vector_argument(y0, "outputs", sys.noutputs)
 
     # Make sure the input arguments match the sizes of the system
     if len(x0) != nstates or \
diff --git a/control/optimal.py b/control/optimal.py
index 31a9998a..3242ac3f 100644
--- a/control/optimal.py
+++ b/control/optimal.py
@@ -35,7 +35,9 @@ import scipy.optimize as opt
 import control as ct
 
 from . import config
+from .config import _process_param, _process_kwargs
 from .iosys import _process_control_disturbance_indices, _process_labels
+from .timeresp import _timeresp_aliases
 
 # Define module default parameter values
 _optimal_trajectory_methods = {'shooting', 'collocation'}
@@ -47,6 +49,19 @@ _optimal_defaults = {
     'optimal.solve_ivp_options': {},
 }
 
+# Parameter and keyword aliases
+_optimal_aliases = {
+    # param:                  ([alias, ...],                [legacy, ...])
+    'integral_cost':          (['trajectory_cost', 'cost'], []),
+    'initial_state':          (['x0', 'X0'],                []),
+    'initial_input':          (['u0', 'U0'],                []),
+    'final_state':            (['xf'],                      []),
+    'final_input':            (['uf'],                      []),
+    'initial_time':           (['T0'],                      []),
+    'trajectory_constraints': (['constraints'],             []),
+    'return_states':          (['return_x'],                []),
+}
+
 
 class OptimalControlProblem():
     """Description of a finite horizon, optimal control problem.
@@ -332,7 +347,8 @@ class OptimalControlProblem():
 
             # Integrate the cost
             costs = np.array(costs)
-           # Approximate the integral using trapezoidal rule
+
+            # Approximate the integral using trapezoidal rule
             cost = np.sum(0.5 * (costs[:-1] + costs[1:]) * dt)
 
         else:
@@ -1018,8 +1034,9 @@ class OptimalControlResult(sp.optimize.OptimizeResult):
 
 # Compute the input for a nonlinear, (constrained) optimal control problem
 def solve_optimal_trajectory(
-        sys, timepts, X0, cost, trajectory_constraints=None,
-        terminal_cost=None, terminal_constraints=None, initial_guess=None,
+        sys, timepts, initial_state=None, integral_cost=None,
+        trajectory_constraints=None, terminal_cost=None,
+        terminal_constraints=None, initial_guess=None,
         basis=None, squeeze=None, transpose=None, return_states=True,
         print_summary=True, log=False, **kwargs):
 
@@ -1044,18 +1061,14 @@ def solve_optimal_trajectory(
     ----------
     sys : `InputOutputSystem`
         I/O system for which the optimal input will be computed.
-
     timepts : 1D array_like
         List of times at which the optimal input should be computed.
-
-    X0 : array_like or number, optional
+    initial_state (or X0) : array_like or number, optional
         Initial condition (default = 0).
-
-    cost : callable
+    integral_cost (or cost) : callable
         Function that returns the integral cost (L) given the current state
-        and input.  Called as ``cost(x, u)``.
-
-    trajectory_constraints : list of tuples, optional
+        and input.  Called as ``integral_cost(x, u)``.
+    trajectory_constraints (or constraints) : list of tuples, optional
         List of constraints that should hold at each point in the time
         vector.  Each element of the list should consist of a tuple with
         first element given by `scipy.optimize.LinearConstraint` or
@@ -1072,52 +1085,23 @@ def solve_optimal_trajectory(
           and compared against the upper and lower bounds.
 
         The constraints are applied at each time point along the trajectory.
-
     terminal_cost : callable, optional
         Function that returns the terminal cost (V) given the final state
         and input.  Called as terminal_cost(x, u).  (For compatibility with
         the form of the cost function, u is passed even though it is often
         not part of the terminal cost.)
-
     terminal_constraints : list of tuples, optional
         List of constraints that should hold at the end of the trajectory.
         Same format as `constraints`.
-
     initial_guess : 1D or 2D array_like
         Initial inputs to use as a guess for the optimal input.  The inputs
         should either be a 2D vector of shape (ninputs, len(timepts)) or a
         1D input of shape (ninputs,) that will be broadcast by extension of
         the time axis.
-
     basis : `BasisFamily`, optional
         Use the given set of basis functions for the inputs instead of
         setting the value of the input at each point in the timepts vector.
 
-    trajectory_method : string, optional
-        Method to use for carrying out the optimization. Currently supported
-        methods are 'shooting' and 'collocation' (continuous time only). The
-        default value is 'shooting' for discrete-time systems and
-        'collocation' for continuous-time systems.
-
-    log : bool, optional
-        If True, turn on logging messages (using Python logging module).
-
-    print_summary : bool, optional
-        If True (default), print a short summary of the computation.
-
-    return_states : bool, optional
-        If True (default), return the values of the state at each time.
-
-    squeeze : bool, optional
-        If True and if the system has a single output, return the system
-        output as a 1D array rather than a 2D array.  If False, return the
-        system output as a 2D array even if the system is SISO.  Default
-        value set by `config.defaults['control.squeeze_time_response']`.
-
-    transpose : bool, optional
-        If True, assume that 2D input arrays are transposed from the standard
-        format.  Used to convert MATLAB-style inputs to our format.
-
     Returns
     -------
     res : `OptimalControlResult`
@@ -1136,8 +1120,27 @@ def solve_optimal_trajectory(
 
     Other Parameters
     ----------------
+    log : bool, optional
+        If True, turn on logging messages (using Python logging module).
     minimize_method : str, optional
         Set the method used by `scipy.optimize.minimize`.
+    print_summary : bool, optional
+        If True (default), print a short summary of the computation.
+    return_states : bool, optional
+        If True (default), return the values of the state at each time.
+    squeeze : bool, optional
+        If True and if the system has a single output, return the system
+        output as a 1D array rather than a 2D array.  If False, return the
+        system output as a 2D array even if the system is SISO.  Default
+        value set by `config.defaults['control.squeeze_time_response']`.
+    trajectory_method : string, optional
+        Method to use for carrying out the optimization. Currently supported
+        methods are 'shooting' and 'collocation' (continuous time only). The
+        default value is 'shooting' for discrete-time systems and
+        'collocation' for continuous-time systems.
+    transpose : bool, optional
+        If True, assume that 2D input arrays are transposed from the standard
+        format.  Used to convert MATLAB-style inputs to our format.
 
     Notes
     -----
@@ -1156,16 +1159,19 @@ def solve_optimal_trajectory(
     `OptimalControlProblem` for more information.
 
     """
-    # Process keyword arguments
-    trajectory_constraints = config._process_legacy_keyword(
-        kwargs, 'constraints', 'trajectory_constraints',
-        trajectory_constraints)
-
-    # Allow 'return_x` as a synonym for 'return_states'
-    return_states = ct.config._get_param(
-        'optimal', 'return_x', kwargs, return_states, pop=True)
-
-    # Process (legacy) method keyword
+    # Process parameter and keyword arguments
+    _process_kwargs(kwargs, _optimal_aliases)
+    X0 = _process_param(
+        'initial_state', initial_state, kwargs, _optimal_aliases, sigval=None)
+    cost = _process_param(
+        'integral_cost', integral_cost, kwargs, _optimal_aliases)
+    trajectory_constraints = _process_param(
+        'trajectory_constraints', trajectory_constraints, kwargs,
+        _optimal_aliases)
+    return_states = _process_param(
+        'return_states', return_states, kwargs, _optimal_aliases, sigval=True)
+
+    # Process (legacy) method keyword (could be minimize or trajectory)
     if kwargs.get('method'):
         method = kwargs.pop('method')
         if method not in _optimal_trajectory_methods:
@@ -1177,7 +1183,8 @@ def solve_optimal_trajectory(
             kwargs['minimize_method'] = method
         else:
             if kwargs.get('trajectory_method'):
-                raise ValueError("'trajectory_method' specified more than once")
+                raise ValueError(
+                    "'trajectory_method' specified more than once")
             warnings.warn(
                 "'method' parameter is deprecated; assuming trajectory_method",
                 FutureWarning)
@@ -1197,8 +1204,8 @@ def solve_optimal_trajectory(
 
 # Create a model predictive controller for an optimal control problem
 def create_mpc_iosystem(
-        sys, timepts, cost, constraints=None, terminal_cost=None,
-        terminal_constraints=None, log=False, **kwargs):
+        sys, timepts, integral_cost=None, trajectory_constraints=None,
+        terminal_cost=None, terminal_constraints=None, log=False, **kwargs):
     """Create a model predictive I/O control system.
 
     This function creates an input/output system that implements a model
@@ -1210,26 +1217,20 @@ def create_mpc_iosystem(
     ----------
     sys : `InputOutputSystem`
         I/O system for which the optimal input will be computed.
-
     timepts : 1D array_like
         List of times at which the optimal input should be computed.
-
-    cost : callable
+    integral_cost (or cost) : callable
         Function that returns the integral cost given the current state
-        and input.  Called as cost(x, u).
-
-    constraints : list of tuples, optional
+        and input.  Called as ``integral_cost(x, u)``.
+    trajectory_constraints (or constraints) : list of tuples, optional
         List of constraints that should hold at each point in the time
         vector.  See `solve_optimal_trajectory` for more details.
-
     terminal_cost : callable, optional
         Function that returns the terminal cost given the final state
         and input.  Called as terminal_cost(x, u).
-
     terminal_constraints : list of tuples, optional
         List of constraints that should hold at the end of the trajectory.
         Same format as `constraints`.
-
     **kwargs
         Additional parameters, passed to `scipy.optimize.minimize` and
         `~control.NonlinearIOSystem`.
@@ -1263,6 +1264,14 @@ def create_mpc_iosystem(
     """
     from .iosys import InputOutputSystem
 
+    # Process parameter and keyword arguments
+    _process_kwargs(kwargs, _optimal_aliases)
+    cost = _process_param(
+        'integral_cost', integral_cost, kwargs, _optimal_aliases)
+    constraints = _process_param(
+        'trajectory_constraints', trajectory_constraints, kwargs,
+        _optimal_aliases)
+
     # Grab the keyword arguments known by this function
     iosys_kwargs = {}
     for kw in InputOutputSystem._kwargs_list:
@@ -1722,17 +1731,17 @@ class OptimalEstimationProblem():
     # Optimal estimate computations
     #
     def compute_estimate(
-            self, Y, U, X0=None, initial_guess=None,
-            squeeze=None, print_summary=True):
+            self, outputs=None, inputs=None, initial_state=None,
+            initial_guess=None, squeeze=None, print_summary=True, **kwargs):
         """Compute the optimal input at state x.
 
         Parameters
         ----------
-        Y : 2D array
+        outputs (or Y) : 2D array
             Measured outputs at each time point.
-        U : 2D array
+        inputs (or U) : 2D array
             Applied inputs at each time point.
-        X0 : 1D array
+        initial_state (or X0) : 1D array
             Expected initial value of the state.
         initial_guess : 2-tuple of 2D arrays
             A 2-tuple consisting of the estimated states and disturbance
@@ -1762,6 +1771,16 @@ class OptimalEstimationProblem():
             Estimated measurement noise for the system trajectory.
 
         """
+        # Argument and keyword processing
+        aliases = _timeresp_aliases | _optimal_aliases
+        _process_kwargs(kwargs, aliases)
+        Y = _process_param('outputs', outputs, kwargs, aliases)
+        U = _process_param('inputs', inputs, kwargs, aliases)
+        X0 = _process_param('initial_state', initial_state, kwargs, aliases)
+
+        if kwargs:
+            raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
         # Store the inputs and outputs (for use in _constraint_function)
         self.u = np.atleast_1d(U).reshape(-1, self.timepts.size)
         self.y = np.atleast_1d(Y).reshape(-1, self.timepts.size)
@@ -1801,7 +1820,6 @@ class OptimalEstimationProblem():
         return OptimalEstimationResult(
             self, res, squeeze=squeeze, print_summary=print_summary)
 
-
     #
     # Create an input/output system implementing an moving horizon estimator
     #
@@ -1809,6 +1827,7 @@ class OptimalEstimationProblem():
     # xhat, u, v, y for all previous time points.  When the system update
     # function is called,
     #
+
     def create_mhe_iosystem(
             self, estimate_labels=None, measurement_labels=None,
             control_labels=None, inputs=None, outputs=None, **kwargs):
@@ -1912,7 +1931,7 @@ class OptimalEstimationProblem():
 
             # Compute the new states and disturbances
             est = self.compute_estimate(
-                Y, U, X0=xhat[:, 0], initial_guess=(xhat, V),
+                Y, U, initial_state=xhat[:, 0], initial_guess=(xhat, V),
                 print_summary=False)
 
             # Restack the new state
@@ -2021,8 +2040,8 @@ class OptimalEstimationResult(sp.optimize.OptimizeResult):
 
 # Compute the finite horizon estimate for a nonlinear system
 def solve_optimal_estimate(
-        sys, timepts, Y, U, trajectory_cost, X0=None,
-        trajectory_constraints=None, initial_guess=None,
+        sys, timepts, outputs=None, inputs=None, integral_cost=None,
+        initial_state=None, trajectory_constraints=None, initial_guess=None,
         squeeze=None, print_summary=True, **kwargs):
 
     """Compute the solution to a finite horizon estimation problem.
@@ -2038,12 +2057,14 @@ def solve_optimal_estimate(
         I/O system for which the optimal input will be computed.
     timepts : 1D array_like
         List of times at which the optimal input should be computed.
-    Y, U : 2D array_like
-        Values of the outputs and inputs at each time point.
-    trajectory_cost : callable
+    outputs (or Y) : 2D array_like
+        Values of the outputs at each time point.
+    inputs (or U) : 2D array_like
+        Values of the inputs at each time point.
+    integral_cost (or cost) : callable
         Function that returns the cost given the current state
         and input.  Called as ``cost(y, u, x0)``.
-    X0 : 1D array_like, optional
+    initial_state (or X0) : 1D array_like, optional
         Mean value of the initial condition (defaults to 0).
     trajectory_constraints : list of tuples, optional
         List of constraints that should hold at each point in the time
@@ -2094,6 +2115,15 @@ def solve_optimal_estimate(
     `OptimalControlProblem` for more information.
 
     """
+    aliases = _timeresp_aliases | _optimal_aliases
+    _process_kwargs(kwargs, aliases)
+    Y = _process_param('outputs', outputs, kwargs, aliases)
+    U = _process_param('inputs', inputs, kwargs, aliases)
+    X0 = _process_param(
+        'initial_state', initial_state, kwargs, aliases)
+    trajectory_cost = _process_param(
+        'integral_cost', integral_cost, kwargs, aliases)
+
     # Set up the optimal control problem
     oep = OptimalEstimationProblem(
         sys, timepts, trajectory_cost,
@@ -2101,7 +2131,7 @@ def solve_optimal_estimate(
 
     # Solve for the optimal input from the current state
     return oep.compute_estimate(
-        Y, U, X0=X0, initial_guess=initial_guess,
+        Y, U, initial_state=X0, initial_guess=initial_guess,
         squeeze=squeeze, print_summary=print_summary)
 
 
@@ -2461,6 +2491,7 @@ def output_range_constraint(sys, lb, ub):
     # Return a nonlinear constraint object based on the polynomial
     return (opt.NonlinearConstraint, _evaluate_output_range_constraint, lb, ub)
 
+
 #
 # Create a constraint on the disturbance input
 #
@@ -2505,6 +2536,7 @@ def disturbance_range_constraint(sys, lb, ub):
 # Utility functions
 #
 
+
 #
 # Process trajectory constraints
 #
@@ -2516,6 +2548,7 @@ def disturbance_range_constraint(sys, lb, ub):
 # internal representation (currently a tuple with the constraint type as the
 # first element.
 #
+
 def _process_constraints(clist, name):
     if clist is None:
         clist = []
@@ -2531,7 +2564,7 @@ def _process_constraints(clist, name):
         if isinstance(constraint, tuple):
             # Original style of constraint
             ctype, fun, lb, ub = constraint
-            if not ctype in [opt.LinearConstraint, opt.NonlinearConstraint]:
+            if ctype not in [opt.LinearConstraint, opt.NonlinearConstraint]:
                 raise TypeError(f"unknown {name} constraint type {ctype}")
             constraint_list.append(constraint)
         elif isinstance(constraint, opt.LinearConstraint):
diff --git a/control/phaseplot.py b/control/phaseplot.py
index 5f5016f5..c17f5fba 100644
--- a/control/phaseplot.py
+++ b/control/phaseplot.py
@@ -977,16 +977,17 @@ def _create_trajectory(
     # Compute the forward trajectory
     if dir == 'forward' or dir == 'both':
         fwdresp = input_output_response(
-            sys, timepts, X0=X0, params=params, ignore_errors=True)
+            sys, timepts, initial_state=X0, params=params, ignore_errors=True)
         if not fwdresp.success and not suppress_warnings:
-            warnings.warn(f"{X0=}, {fwdresp.message}")
+            warnings.warn(f"initial_state={X0}, {fwdresp.message}")
 
     # Compute the reverse trajectory
     if dir == 'reverse' or dir == 'both':
         revresp = input_output_response(
-            revsys, timepts, X0=X0, params=params, ignore_errors=True)
+            revsys, timepts, initial_state=X0, params=params,
+            ignore_errors=True)
         if not revresp.success and not suppress_warnings:
-            warnings.warn(f"{X0=}, {revresp.message}")
+            warnings.warn(f"initial_state={X0}, {revresp.message}")
 
     # Create the trace to plot
     if dir == 'forward':
diff --git a/control/statesp.py b/control/statesp.py
index 2eba44df..65529b99 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -1525,7 +1525,7 @@ class StateSpace(NonlinearIOSystem, LTI):
             return (self.C @ x).reshape((-1,)) \
                 + (self.D @ u).reshape((-1,))  # return as row vector
 
-    # convenience alias, import needs to go over the submodule to avoid circular imports
+    # convenience alias, import needs submodule to avoid circular imports
     initial_response = control.timeresp.initial_response
 
 
diff --git a/control/timeresp.py b/control/timeresp.py
index d708dd2e..bd549589 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -45,6 +45,7 @@ from numpy import einsum, maximum, minimum
 from scipy.linalg import eig, eigvals, matrix_balance, norm
 
 from . import config
+from . config import _process_kwargs, _process_param
 from .exception import pandas_check
 from .iosys import NamedSignal, isctime, isdtime
 from .timeplot import time_response_plot
@@ -53,6 +54,21 @@ __all__ = ['forced_response', 'step_response', 'step_info',
            'initial_response', 'impulse_response', 'TimeResponseData',
            'TimeResponseList']
 
+# Dictionary of aliases for time response commands
+_timeresp_aliases = {
+    # param:            ([alias, ...], [legacy, ...])
+    'timepts':          (['T'],        []),
+    'inputs':           (['U'],        ['u']),
+    'outputs':          (['Y'],        ['y']),
+    'initial_state':    (['X0'],       ['x0']),
+    'final_output':     (['yfinal'],   []),
+    'return_states':    (['return_x'], []),
+    'evaluation_times': (['t_eval'],   []),
+    'timepts_num':      (['T_num'],    []),
+    'input_indices':    (['input'],    []),
+    'output_indices':   (['output'],   []),
+}
+
 
 class TimeResponseData:
     """Input/output system time response data.
@@ -197,7 +213,7 @@ class TimeResponseData:
     names of the appropriate signals::
 
       sys = ct.rss(4, 2, 1)
-      resp = ct.initial_response(sys, X0=[1, 1, 1, 1])
+      resp = ct.initial_response(sys, initial_state=[1, 1, 1, 1])
       plt.plot(resp.time, resp.outputs['y[0]'])
 
     In the case of multi-trace data, the responses should be indexed using
@@ -343,7 +359,7 @@ class TimeResponseData:
             # Make sure the shape is OK
             if multi_trace and \
                (self.x.ndim != 3 or self.x.shape[1] != self.ntraces) or \
-               not multi_trace and self.x.ndim != 2 :
+               not multi_trace and self.x.ndim != 2:
                 raise ValueError("State vector is the wrong shape")
 
             # Make sure time dimension of state is the right length
@@ -700,6 +716,7 @@ class TimeResponseData:
         """
         return time_response_plot(self, *args, **kwargs)
 
+
 #
 # Time response data list class
 #
@@ -900,8 +917,10 @@ def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
 
 
 # Forced response of a linear system
-def forced_response(sysdata, T=None, U=0., X0=0., transpose=False, params=None,
-                    interpolate=False, return_x=None, squeeze=None):
+def forced_response(
+        sysdata, timepts=None, inputs=0., initial_state=0., transpose=False,
+        params=None, interpolate=False, return_states=None, squeeze=None,
+        **kwargs):
     """Compute the output of a linear system given the input.
 
     As a convenience for parameters `U`, `X0`: Numbers (scalars) are
@@ -915,46 +934,36 @@ def forced_response(sysdata, T=None, U=0., X0=0., transpose=False, params=None,
     ----------
     sysdata : I/O system or list of I/O systems
         I/O system(s) for which forced response is computed.
-
-    T : array_like, optional for discrete LTI `sys`
+    timepts (or T) : array_like, optional for discrete LTI `sys`
         Time steps at which the input is defined; values must be evenly
-        spaced.  If None, `U` must be given and ``len(U)`` time steps of
-        sys.dt are simulated. If sys.dt is None or True (undetermined
-        time step), a time step of 1.0 is assumed.
-
-    U : array_like or float, optional
-        Input array giving input at each time `T`.  If `U` is None or 0,
-        `T` must be given, even for discrete-time systems. In this case,
-        for continuous-time systems, a direct calculation of the matrix
-        exponential is used, which is faster than the general interpolating
-        algorithm used otherwise.
-
-    X0 : array_like or float, default=0.
+        spaced.  If None, `inputs` must be given and ``len(inputs)`` time
+        steps of `sys.dt` are simulated. If `sys.dt` is None or True
+        (undetermined time step), a time step of 1.0 is assumed.
+    inputs (or U) : array_like or float, optional
+        Input array giving input at each time in `timepts`.  If `inputs` is
+        None or 0, `timepts` must be given, even for discrete-time
+        systems. In this case, for continuous-time systems, a direct
+        calculation of the matrix exponential is used, which is faster than
+        the general interpolating algorithm used otherwise.
+    initial_state (or X0) : array_like or float, default=0.
         Initial condition.
-
     params : dict, optional
         If system is a nonlinear I/O system, set parameter values.
-
     transpose : bool, default=False
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and `scipy.signal.lsim`).
-
     interpolate : bool, default=False
         If True and system is a discrete-time system, the input will
         be interpolated between the given time steps and the output
         will be given at system sampling rate.  Otherwise, only return
         the output at the times given in `T`.  No effect on continuous
         time simulations.
-
-    return_x : bool, default=None
-        Used if the time response data is assigned to a tuple:
-
-            * If False, return only the time and output vectors.
-            * If True, also return the the state vector.
-            * If None, determine the returned variables by
-              `config.defaults['forced_response.return_x']`, which was True
-              before version 0.9 and is False since then.
-
+    return_states (or return_x) : bool, default=None
+        Used if the time response data is assigned to a tuple.  If False,
+        return only the time and output vectors.  If True, also return the
+        the state vector.  If None, determine the returned variables by
+        `config.defaults['forced_response.return_x']`, which was True
+        before version 0.9 and is False since then.
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then
         the output response is returned as a 1D array (indexed by time).
@@ -984,8 +993,8 @@ def forced_response(sysdata, T=None, U=0., X0=0., transpose=False, params=None,
         True, the array is 1D (indexed by time).  If the system is not SISO or
         `squeeze` is False, the array is 2D (indexed by output and time).
     resp.states : array
-        Time evolution of the state vector, represented as a 2D array indexed by
-        state and time.
+        Time evolution of the state vector, represented as a 2D array
+        indexed by state and time.
     resp.inputs : array
         Input(s) to the system, indexed by input and time.
 
@@ -1016,8 +1025,9 @@ def forced_response(sysdata, T=None, U=0., X0=0., transpose=False, params=None,
     Examples
     --------
     >>> G = ct.rss(4)
-    >>> T = np.linspace(0, 10)
-    >>> T, yout = ct.forced_response(G, T=T)
+    >>> timepts = np.linspace(0, 10)
+    >>> inputs = np.sin(timepts)
+    >>> tout, yout = ct.forced_response(G, timepts, inputs)
 
     See :ref:`time-series-convention` and
     :ref:`package-configuration-parameters`.
@@ -1027,13 +1037,26 @@ def forced_response(sysdata, T=None, U=0., X0=0., transpose=False, params=None,
     from .statesp import StateSpace, _convert_to_statespace
     from .xferfcn import TransferFunction
 
+    # Process keyword arguments
+    _process_kwargs(kwargs, _timeresp_aliases)
+    T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+    U = _process_param('inputs', inputs, kwargs, _timeresp_aliases, sigval=0.)
+    X0 = _process_param(
+        'initial_state', initial_state, kwargs, _timeresp_aliases, sigval=0.)
+    return_x = _process_param(
+        'return_states', return_states, kwargs, _timeresp_aliases, sigval=None)
+
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
     # If passed a list, recursively call individual responses with given T
     if isinstance(sysdata, (list, tuple)):
         responses = []
         for sys in sysdata:
             responses.append(forced_response(
-                sys, T, U=U, X0=X0, transpose=transpose, params=params,
-                interpolate=interpolate, return_x=return_x, squeeze=squeeze))
+                sys, T, inputs=U, initial_state=X0, transpose=transpose,
+                params=params, interpolate=interpolate,
+                return_states=return_x, squeeze=squeeze))
         return TimeResponseList(responses)
     else:
         sys = sysdata
@@ -1309,8 +1332,9 @@ def _process_time_response(
 
 
 def step_response(
-        sysdata, T=None, X0=0, input=None, output=None, T_num=None,
-        transpose=False, return_x=False, squeeze=None, params=None):
+        sysdata, timepts=None, initial_state=0., input_indices=None,
+        output_indices=None, timepts_num=None, transpose=False,
+        return_states=False, squeeze=None, params=None, **kwargs):
     # pylint: disable=W0622
     """Compute the step response for a linear system.
 
@@ -1327,8 +1351,7 @@ def step_response(
     ----------
     sysdata : I/O system or list of I/O systems
         I/O system(s) for which step response is computed.
-
-    T : array_like or float, optional
+    timepts (or T) : array_like or float, optional
         Time vector, or simulation time duration if a number. If `T` is not
         provided, an attempt is made to create it automatically from the
         dynamics of the system. If the system continuous time, the time
@@ -1338,37 +1361,29 @@ def step_response(
         this results in too many time steps (>5000), dt is reduced. If the
         system is discrete time, only tfinal is computed, and final is
         reduced if it requires too many simulation steps.
-
-    X0 : array_like or float, optional
+    initial_state (or X0) : array_like or float, optional
         Initial condition (default = 0).  This can be used for a nonlinear
         system where the origin is not an equilibrium point.
-
-    input : int, optional
+    input_indices (or input) : int or list of int, optional
         Only compute the step response for the listed input.  If not
         specified, the step responses for each independent input are
         computed (as separate traces).
-
-    output : int, optional
+    output_indices (or output) : int, optional
         Only report the step response for the listed output.  If not
         specified, all outputs are reported.
-
     params : dict, optional
         If system is a nonlinear I/O system, set parameter values.
-
-    T_num : int, optional
+    timepts_num (or T_num) : int, optional
         Number of time steps to use in simulation if `T` is not provided as
         an array (auto-computed if not given); ignored if the system is
         discrete time.
-
     transpose : bool, optional
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and `scipy.signal.lsim`).  Default
         value is False.
-
-    return_x : bool, optional
+    return_states (or return_x) : bool, optional
         If True, return the state vector when assigning to a tuple
         (default = False).  See `forced_response` for more details.
-
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then
         the output response is returned as a 1D array (indexed by time).
@@ -1405,6 +1420,24 @@ def step_response(
     from .statesp import _convert_to_statespace
     from .xferfcn import TransferFunction
 
+    # Process keyword arguments
+    _process_kwargs(kwargs, _timeresp_aliases)
+    T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+    X0 = _process_param(
+        'initial_state', initial_state, kwargs, _timeresp_aliases, sigval=0.)
+    input = _process_param(
+        'input_indices', input_indices, kwargs, _timeresp_aliases)
+    output = _process_param(
+        'output_indices', output_indices, kwargs, _timeresp_aliases)
+    return_x = _process_param(
+        'return_states', return_states, kwargs, _timeresp_aliases,
+        sigval=False)
+    T_num = _process_param(
+        'timepts_num', timepts_num, kwargs, _timeresp_aliases)
+
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
     # Create the time and input vectors
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=True)
@@ -1417,8 +1450,9 @@ def step_response(
         responses = []
         for sys in sysdata:
             responses.append(step_response(
-                sys, T, X0=X0, input=input, output=output, T_num=T_num,
-                transpose=transpose, return_x=return_x, squeeze=squeeze,
+                sys, T, initial_state=X0, input_indices=input,
+                output_indices=output, timepts_num=T_num,
+                transpose=transpose, return_states=return_x, squeeze=squeeze,
                 params=params))
         return TimeResponseList(responses)
     else:
@@ -1435,6 +1469,21 @@ def step_response(
     if isinstance(sys, LTI) and sys.nstates is None:
         sys = _convert_to_statespace(sys)
 
+    # Only single input and output are allowed for now
+    if isinstance(input, (list, tuple)):
+        if len(input_indices) > 1:
+            raise NotImplementedError("list of input indices not allowed")
+        input = input[0]
+    elif isinstance(input, str):
+        raise NotImplementedError("named inputs not allowed")
+
+    if isinstance(output, (list, tuple)):
+        if len(output_indices) > 1:
+            raise NotImplementedError("list of output indices not allowed")
+        output = output[0]
+    elif isinstance(output, str):
+        raise NotImplementedError("named outputs not allowed")
+
     # Set up arrays to handle the output
     ninputs = sys.ninputs if input is None else 1
     noutputs = sys.noutputs if output is None else 1
@@ -1482,8 +1531,10 @@ def step_response(
         trace_types=trace_types, plot_inputs=False)
 
 
-def step_info(sysdata, T=None, T_num=None, yfinal=None, params=None,
-              SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)):
+def step_info(
+        sysdata, timepts=None, timepts_num=None, final_output=None,
+        params=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9),
+        **kwargs):
     """Step response characteristics (rise time, settling time, etc).
 
     Parameters
@@ -1491,15 +1542,15 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None, params=None,
     sysdata : `StateSpace` or `TransferFunction` or array_like
         The system data. Either LTI system to simulate (`StateSpace`,
         `TransferFunction`), or a time series of step response data.
-    T : array_like or float, optional
+    timepts (or T) : array_like or float, optional
         Time vector, or simulation time duration if a number (time vector is
         auto-computed if not given, see `step_response` for more detail).
         Required, if sysdata is a time series of response data.
-    T_num : int, optional
+    timepts_num (or T_num) : int, optional
         Number of time steps to use in simulation if `T` is not provided as
         an array; auto-computed if not given; ignored if sysdata is a
         discrete-time system or a time series or response data.
-    yfinal : scalar or array_like, optional
+    final_output (or yfinal) : scalar or array_like, optional
         Steady-state response. If not given, sysdata.dcgain() is used for
         systems to simulate and the last value of the the response data is
         used for a given time series of response data. Scalar for SISO,
@@ -1581,8 +1632,20 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None, params=None,
     from .statesp import StateSpace
     from .xferfcn import TransferFunction
 
+    # Process keyword arguments
+    _process_kwargs(kwargs, _timeresp_aliases)
+    T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+    T_num = _process_param(
+        'timepts_num', timepts_num, kwargs, _timeresp_aliases)
+    yfinal = _process_param(
+        'final_output', final_output, kwargs, _timeresp_aliases)
+
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
     if isinstance(sysdata, (StateSpace, TransferFunction, NonlinearIOSystem)):
-        T, Yout = step_response(sysdata, T, squeeze=False, params=params)
+        T, Yout = step_response(
+            sysdata, T, timepts_num=T_num, squeeze=False, params=params)
         if yfinal:
             InfValues = np.atleast_2d(yfinal)
         else:
@@ -1700,8 +1763,9 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None, params=None,
 
 
 def initial_response(
-        sysdata, T=None, X0=0, output=None, T_num=None, params=None,
-        transpose=False, return_x=False, squeeze=None):
+        sysdata, timepts=None, initial_state=0, output_indices=None,
+        timepts_num=None, params=None, transpose=False, return_states=False,
+        squeeze=None, **kwargs):
     # pylint: disable=W0622
     """Compute the initial condition response for a linear system.
 
@@ -1716,39 +1780,28 @@ def initial_response(
     ----------
     sysdata : I/O system or list of I/O systems
         I/O system(s) for which initial response is computed.
-
-    sys : `StateSpace` or `TransferFunction`
-        LTI system to simulate.
-
-    T :  array_like or float, optional
+    timepts (or T) :  array_like or float, optional
         Time vector, or simulation time duration if a number (time vector is
         auto-computed if not given; see  `step_response` for more detail).
-
-    X0 : array_like or float, optional
+    initial_state (or X0) : array_like or float, optional
         Initial condition (default = 0).  Numbers are converted to constant
         arrays with the correct shape.
-
-    output : int
+    output_indices (or output) : int
         Index of the output that will be used in this simulation. Set
         to None to not trim outputs.
-
-    T_num : int, optional
-        Number of time steps to use in simulation if `T` is not provided as
-        an array (auto-computed if not given); ignored if the system is
-        discrete time.
-
+    timepts_num (or T_num) : int, optional
+        Number of time steps to use in simulation if `timepts` is not
+        provided as an array (auto-computed if not given); ignored if the
+        system is discrete time.
     params : dict, optional
         If system is a nonlinear I/O system, set parameter values.
-
     transpose : bool, optional
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and `scipy.signal.lsim`).  Default
         value is False.
-
-    return_x : bool, optional
+    return_states (or return_x) : bool, optional
         If True, return the state vector when assigning to a tuple
         (default = False).  See `forced_response` for more details.
-
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then
         the output response is returned as a 1D array (indexed by time).
@@ -1781,6 +1834,22 @@ def initial_response(
     >>> T, yout = ct.initial_response(G)
 
     """
+    # Process keyword arguments
+    _process_kwargs(kwargs, _timeresp_aliases)
+    T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+    X0 = _process_param(
+        'initial_state', initial_state, kwargs, _timeresp_aliases, sigval=0.)
+    output = _process_param(
+        'output_indices', output_indices, kwargs, _timeresp_aliases)
+    return_x = _process_param(
+        'return_states', return_states, kwargs, _timeresp_aliases,
+        sigval=False)
+    T_num = _process_param(
+        'timepts_num', timepts_num, kwargs, _timeresp_aliases)
+
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
     # Create the time and input vectors
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=False)
@@ -1793,8 +1862,9 @@ def initial_response(
         responses = []
         for sys in sysdata:
             responses.append(initial_response(
-                sys, T, X0=X0, output=output, T_num=T_num, transpose=transpose,
-                return_x=return_x, squeeze=squeeze, params=params))
+                sys, T, initial_state=X0, output_indices=output,
+                timepts_num=T_num, transpose=transpose,
+                return_states=return_x, squeeze=squeeze, params=params))
         return TimeResponseList(responses)
     else:
         sys = sysdata
@@ -1820,8 +1890,9 @@ def initial_response(
 
 
 def impulse_response(
-        sysdata, T=None, input=None, output=None, T_num=None,
-        transpose=False, return_x=False, squeeze=None):
+        sysdata, timepts=None, input_indices=None, output_indices=None,
+        timepts_num=None, transpose=False, return_states=False, squeeze=None,
+        **kwargs):
     # pylint: disable=W0622
     """Compute the impulse response for a linear system.
 
@@ -1838,34 +1909,27 @@ def impulse_response(
     ----------
     sysdata : I/O system or list of I/O systems
         I/O system(s) for which impulse response is computed.
-
-    T : array_like or float, optional
+    timepts (or T) : array_like or float, optional
         Time vector, or simulation time duration if a scalar (time vector is
         auto-computed if not given; see `step_response` for more detail).
-
-    input : int, optional
+    input_indices (or input) : int, optional
         Only compute the impulse response for the listed input.  If not
         specified, the impulse responses for each independent input are
         computed.
-
-    output : int, optional
+    output_indices (or output) : int, optional
         Only report the step response for the listed output.  If not
         specified, all outputs are reported.
-
-    T_num : int, optional
+    timepts_num (or T_num) : int, optional
         Number of time steps to use in simulation if `T` is not provided as
         an array (auto-computed if not given); ignored if the system is
         discrete time.
-
     transpose : bool, optional
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and `scipy.signal.lsim`).  Default
         value is False.
-
-    return_x : bool, optional
+    return_states (or return_x) : bool, optional
         If True, return the state vector when assigning to a tuple
         (default = False).  See `forced_response` for more details.
-
     squeeze : bool, optional
         By default, if a system is single-input, single-output (SISO) then
         the output response is returned as a 1D array (indexed by time).
@@ -1904,6 +1968,22 @@ def impulse_response(
     from .lti import LTI
     from .statesp import _convert_to_statespace
 
+    # Process keyword arguments
+    _process_kwargs(kwargs, _timeresp_aliases)
+    T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+    input = _process_param(
+        'input_indices', input_indices, kwargs, _timeresp_aliases)
+    output = _process_param(
+        'output_indices', output_indices, kwargs, _timeresp_aliases)
+    return_x = _process_param(
+        'return_states', return_states, kwargs, _timeresp_aliases,
+        sigval=False)
+    T_num = _process_param(
+        'timepts_num', timepts_num, kwargs, _timeresp_aliases)
+
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): ", str(kwargs))
+
     # Create the time and input vectors
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=False)
@@ -1937,6 +2017,20 @@ def impulse_response(
                       "output.\n"
                       "Results may be meaningless!")
 
+    # Only single input and output are allowed for now
+    if isinstance(input, (list, tuple)):
+        if len(input_indices) > 1:
+            raise NotImplementedError("list of input indices not allowed")
+        input = input[0]
+    elif isinstance(input, str):
+        raise NotImplementedError("named inputs not allowed")
+
+    if isinstance(output, (list, tuple)):
+        if len(output_indices) > 1:
+            raise NotImplementedError("list of output indices not allowed")
+        output = output[0]
+    elif isinstance(output, str):
+        raise NotImplementedError("named outputs not allowed")
 
     # Set up arrays to handle the output
     ninputs = sys.ninputs if input is None else 1
diff --git a/doc/develop.rst b/doc/develop.rst
index 6386d605..c9b6738a 100644
--- a/doc/develop.rst
+++ b/doc/develop.rst
@@ -217,7 +217,7 @@ frequency, etc:
 
     def model_update(t, x, u, params)
     resp = initial_response(sys, timepts, x0)             # x0 required
-    resp = input_output_response(sys, timepts, u, x0=x0)  # u required
+    resp = input_output_response(sys, timepts, u, x0)     # u required
     resp = TimeResponseData(
         timepts, outputs, states=states, inputs=inputs)
 
@@ -301,6 +301,157 @@ Time and frequency responses:
       mag, phase, omega = fresp.magnitude, fresp.phase, fresp.omega
 
 
+Parameter aliases
+-----------------
+
+As described above, parameter names are generally longer strings that
+describe the purpose of the parameter.  Similar to `matplotlib` (e.g.,
+the use of `lw` as an alias for `linewidth`), some commonly used
+parameter names can be specified using an "alias" that allows the use
+of a shorter key.
+
+Named parameter and keyword variable aliases are processed using the
+:func:`config._process_kwargs` and :func:`config._process_param`
+functions.  These functions allow the specification of a list of
+aliases and a list of legacy keys for a given named parameter or
+keyword.  To make use of these functions, the
+:func:`~config._process_kwargs` is first called to update the `kwargs`
+variable by replacing aliases with the full key::
+
+  _process_kwargs(kwargs, aliases)
+
+The values for named parameters can then be assigned to a local
+variable using a call to :func:`~config._process_param` of the form::
+
+  var = _process_param('param', param, kwargs, aliases)
+
+where `param` is the named parameter used in the function signature
+and var is the local variable in the function (may also be `param`,
+but doesn't have to be).
+
+For example, the following structure is used in `input_output_response`::
+
+  def input_output_response(
+          sys, timepts=None, inputs=0., initial_state=0., params=None,
+          ignore_errors=False, transpose=False, return_states=False,
+          squeeze=None, solve_ivp_kwargs=None, evaluation_times='T', **kwargs):
+      """Compute the output response of a system to a given input.
+
+      ... rest of docstring ...
+
+      """
+      _process_kwargs(kwargs, _timeresp_aliases)
+      T = _process_param('timepts', timepts, kwargs, _timeresp_aliases)
+      U = _process_param('inputs', inputs, kwargs, _timeresp_aliases, sigval=0.)
+      X0 = _process_param(
+          'initial_state', initial_state, kwargs, _timeresp_aliases, sigval=0.)
+
+Note that named parameters that have a default value other than None
+must given the signature value (`sigval`) so that
+`~config._process_param` can detect if the value has been set (and
+issue an error if there is an attempt to set the value multiple times
+using alias or legacy keys).
+
+The alias mapping is a dictionary that returns a tuple consisting of
+valid aliases and legacy aliases::
+
+  alias_mapping = {
+      'argument_name_1': (['alias', ...], ['legacy', ...]),
+       ...}
+
+If an alias is present in the dictionary of keywords, it will be used
+to set the value of the argument.  If a legacy keyword is used, a
+warning is issued.
+
+The following tables summarize the aliases that are currently in use
+through the python-control package:
+
+Time response aliases (via `timeresp._timeresp_aliases`):
+
+  .. list-table::
+     :header-rows: 1
+
+     * - Key
+       - Aliases
+       - Legacy keys
+       - Comment
+     * - evaluation_times
+       - t_eval
+       -
+       - List of times to evaluate the time response (defaults to `timepts`).
+     * - final_output
+       - yfinal
+       -
+       - Final value of the output (used for :func:`step_info`)
+     * - initial_state
+       - X0
+       - x0
+       - Initial value of the state variable.
+     * - input_indices
+       - input
+       -
+       - Index(es) to use for the input (used in
+         :func:`step_response`, :func:`impulse_response`.
+     * - inputs
+       - U
+       - u
+       - Value(s) of the input variable (time trace or individual point).
+     * - output_indices
+       - output
+       -
+       - Index(es) to use for the output (used in
+         :func:`step_response`, :func:`impulse_response`.
+     * - outputs
+       - Y
+       - y
+       - Value(s) of the output variable (time trace or individual point).
+     * - return_states
+       - return_x
+       -
+       - Return the state when accessing a response via a tuple.
+     * - timepts
+       - T
+       -
+       - List of time points for time response functions.
+     * - timepts_num
+       - T_num
+       -
+       - Number of points to use (e.g., if `timepts` is just the final time).
+
+Optimal control aliases (via `optimal._optimal_aliases`:
+
+  .. list-table::
+     :header-rows: 1
+
+     * - Key
+       - Aliases
+       - Comment
+     * - final_state
+       - xf
+       - Final state for trajectory generation problems (flatsys, optimal).
+     * - final_input
+       - uf
+       - Final input for trajectory generation problems (flatsys).
+     * - initial_state
+       - x0, X0
+       - Initial state for optimization problems (flatsys, optimal).
+     * - initial_input
+       - u0, U0
+       - Initial input for trajectory generation problems (flatsys).
+     * - initial_time
+       - T0
+       - Initial time for optimization problems.
+     * - integral_cost
+       - trajectory_cost, cost
+       - Cost function that is integrated along a trajectory.
+     * - return_states
+       - return_x
+       - Return the state when accessing a response via a tuple.
+     * - trajectory_constraints
+       - constraints
+       - List of constraints that hold along a trajectory (flatsys, optimal)
+
+
 Documentation Guidelines
 ========================
 
@@ -625,6 +776,8 @@ processing and parsing operations:
    :toctree: generated/
 
    config._process_legacy_keyword
+   config._process_kwargs
+   config._process_param
    exception.cvxopt_check
    exception.pandas_check
    exception.slycot_check
diff --git a/doc/iosys.rst b/doc/iosys.rst
index c1d58fb6..5e51e7f0 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -269,7 +269,8 @@ Finally, we simulate the closed loop system:
 
   # Simulate the system
   Ld = 30
-  resp = ct.input_output_response(closed, timepts, U=Ld, X0=[15, 20])
+  resp = ct.input_output_response(
+      closed, timepts, inputs=Ld, initial_state=[15, 20])
   cplt = resp.plot(
       plot_inputs=False, overlay_signals=True, legend_loc='upper left')
   cplt.axes[0, 0].axhline(Ld, linestyle='--', color='black')
diff --git a/doc/optimal.rst b/doc/optimal.rst
index 07131b53..41625689 100644
--- a/doc/optimal.rst
+++ b/doc/optimal.rst
@@ -294,7 +294,7 @@ Plotting the results:
   # Simulate the system dynamics (open loop)
   resp = ct.input_output_response(
       vehicle, timepts, result.inputs, x0,
-      t_eval=np.linspace(0, Tf, 100))
+      evaluation_times=np.linspace(0, Tf, 100))
   t, y, u = resp.time, resp.outputs, resp.inputs
 
   plt.subplot(3, 1, 1)
diff --git a/doc/stochastic.rst b/doc/stochastic.rst
index b09213f1..881cf234 100644
--- a/doc/stochastic.rst
+++ b/doc/stochastic.rst
@@ -375,7 +375,8 @@ Given noisy measurements :math:`y` and control inputs :math:`u`, an
 estimate of the states over the time points can be computed using the
 :func:`~optimal.OptimalEstimationProblem.compute_estimate` method::
 
-  estim = oep.compute_optimal(Y, U[, X0=x0, initial_guess=(xhat, v)])
+  estim = oep.compute_optimal(
+      Y, U[, initial_state=x0, initial_guess=(xhat, v)])
   xhat, v, w = estim.states, estim.inputs, estim.outputs
 
 For discrete-time systems, the


Test patch:
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 3b7de8b8..aa243b60 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -46,7 +46,7 @@ function_skiplist = [
 
 # List of keywords that we can skip testing (special cases)
 keyword_skiplist = {
-    control.input_output_response: ['method'],
+    control.input_output_response: ['method', 't_eval'],    # solve_ivp_kwargs
     control.nyquist_plot: ['color'],                        # separate check
     control.optimal.solve_optimal_trajectory:
       ['method', 'return_x'],                               # deprecated
@@ -637,14 +637,15 @@ def _check_parameter_docs(
         docstring = docstring[start:]
 
     # Look for the parameter name in the docstring
+    argname_ = argname + r"( \(or .*\))*"
     if match := re.search(
-            "\n" + r"((\w+|\.{3}), )*" + argname + r"(, (\w+|\.{3}))*:",
+            "\n" + r"((\w+|\.{3}), )*" + argname_ + r"(, (\w+|\.{3}))*:",
             docstring):
         # Found the string, but not in numpydoc form
         _warn(f"{funcname}: {argname} docstring missing space")
 
     elif not (match := re.search(
-            "\n" + r"((\w+|\.{3}), )*" + argname + r"(, (\w+|\.{3}))* :",
+            "\n" + r"((\w+|\.{3}), )*" + argname_ + r"(, (\w+|\.{3}))* :",
             docstring)):
         if fail_if_missing:
             _fail(f"{funcname} '{argname}' not documented")
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 1eb1a1fd..10eb7fb6 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2161,7 +2161,8 @@ def test_operating_point():
     assert isinstance(op_point[1], np.ndarray)
     assert isinstance(op_point[2], np.ndarray)
 
-    with pytest.warns(FutureWarning, match="return_outputs"):
+    with pytest.warns(
+            (FutureWarning, PendingDeprecationWarning), match="return_outputs"):
         op_point = ct.find_operating_point(sys, 0, 0, return_y=True)
         assert len(op_point) == 3
         assert isinstance(op_point[0], np.ndarray)
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 2e491900..54a1fc76 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -30,6 +30,7 @@ import control.tests.statesp_test as statesp_test
 import control.tests.statefbk_test as statefbk_test
 import control.tests.stochsys_test as stochsys_test
 import control.tests.timeplot_test as timeplot_test
+import control.tests.timeresp_test as timeresp_test
 import control.tests.trdata_test as trdata_test
 
 
@@ -260,9 +261,12 @@ kwarg_unittest = {
     'find_eqpt': iosys_test.test_find_operating_point,
     'find_operating_point': iosys_test.test_find_operating_point,
     'flatsys.flatsys': test_unrecognized_kwargs,
+    'forced_response': timeresp_test.test_timeresp_aliases,
     'frd': frd_test.TestFRD.test_unrecognized_keyword,
     'gangof4': test_matplotlib_kwargs,
     'gangof4_plot': test_matplotlib_kwargs,
+    'impulse_response': timeresp_test.test_timeresp_aliases,
+    'initial_response': timeresp_test.test_timeresp_aliases,
     'input_output_response': test_unrecognized_kwargs,
     'interconnect': interconnect_test.test_interconnect_exceptions,
     'time_response_plot': timeplot_test.test_errors,
@@ -294,6 +298,8 @@ kwarg_unittest = {
     'set_defaults': test_unrecognized_kwargs,
     'singular_values_plot': test_matplotlib_kwargs,
     'ss': test_unrecognized_kwargs,
+    'step_info': timeresp_test.test_timeresp_aliases,
+    'step_response': timeresp_test.test_timeresp_aliases,
     'LTI.to_ss': test_unrecognized_kwargs,          # tested via 'ss'
     'ss2io': test_unrecognized_kwargs,
     'ss2tf': test_unrecognized_kwargs,
@@ -342,6 +348,7 @@ kwarg_unittest = {
     'NonlinearIOSystem.__init__':
         interconnect_test.test_interconnect_exceptions,
     'StateSpace.__init__': test_unrecognized_kwargs,
+    'StateSpace.initial_response': timeresp_test.test_timeresp_aliases,
     'StateSpace.sample': test_unrecognized_kwargs,
     'TimeResponseData.__call__': trdata_test.test_response_copy,
     'TimeResponseData.plot': timeplot_test.test_errors,
@@ -356,6 +363,8 @@ kwarg_unittest = {
         optimal_test.test_ocp_argument_errors,
     'optimal.OptimalEstimationProblem.__init__':
         optimal_test.test_oep_argument_errors,
+    'optimal.OptimalEstimationProblem.compute_estimate':
+        stochsys_test.test_oep,
     'optimal.OptimalEstimationProblem.create_mhe_iosystem':
         optimal_test.test_oep_argument_errors,
     'phaseplot.streamlines': test_matplotlib_kwargs,
diff --git a/control/tests/phaseplot_test.py b/control/tests/phaseplot_test.py
index 106dee6f..6ea6411d 100644
--- a/control/tests/phaseplot_test.py
+++ b/control/tests/phaseplot_test.py
@@ -167,7 +167,8 @@ def test_phaseplane_errors():
         lambda t, x, u, params: np.array(
             [x[1], -0.25 * (x[0] - 0.01 * x[0]**3) - 0.1 * x[1]]),
         states=2, inputs=0)
-    with pytest.warns(UserWarning, match=r"X0=array\(.*\), solve_ivp failed"):
+    with pytest.warns(
+            UserWarning, match=r"initial_state=\[.*\], solve_ivp failed"):
         ct.phase_plane_plot(
             sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
             plot_separatrices=False)
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index b0ddea61..2f3aa512 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -1638,9 +1638,9 @@ def test_convenience_aliases():
     # Make sure that unrecognized keywords for response functions are caught
     for method in [LTI.impulse_response, LTI.initial_response,
                    LTI.step_response]:
-        with pytest.raises(TypeError, match="unexpected keyword"):
+        with pytest.raises(TypeError, match="unrecognized keyword"):
             method(sys, unknown=True)
-    with pytest.raises(TypeError, match="unexpected keyword"):
+    with pytest.raises(TypeError, match="unrecognized keyword"):
         LTI.forced_response(sys, [0, 1], [1, 1], unknown=True)
 
 
diff --git a/control/tests/stochsys_test.py b/control/tests/stochsys_test.py
index 0bbf49b5..4c0d9665 100644
--- a/control/tests/stochsys_test.py
+++ b/control/tests/stochsys_test.py
@@ -404,6 +404,10 @@ def test_oep(dt):
     np.testing.assert_allclose(
         est3.states[:, -1], res3.states[:, -1], atol=meas_mag, rtol=meas_mag)
 
+    # Make sure unknown keywords generate an error
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        est3 = oep1.compute_estimate(Y3, U, unknown=True)
+
 
 @pytest.mark.slow
 def test_mhe():
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index a410bf30..aa498720 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -341,10 +341,15 @@ class TestTimeresp:
         t = tsystem.t
         yref = tsystem.ystep
         _t, y_00 = step_response(sys, T=t, input=0, output=0)
-        _t, y_11 = step_response(sys, T=t, input=1, output=1)
         np.testing.assert_array_almost_equal(y_00, yref, decimal=4)
+
+        _t, y_11 = step_response(sys, T=t, input=1, output=1)
         np.testing.assert_array_almost_equal(y_11, yref, decimal=4)
 
+        _t, y_01 = step_response(
+            sys, T=t, input_indices=[0], output_indices=[1])
+        np.testing.assert_array_almost_equal(y_01, 0 * yref, decimal=4)
+
         # Make sure we get the same result using MIMO step response
         response = step_response(sys, T=t)
         np.testing.assert_allclose(response.y[0, 0, :], y_00)
@@ -354,6 +359,11 @@ class TestTimeresp:
         np.testing.assert_allclose(response.u[0, 1, :], 0)
         np.testing.assert_allclose(response.u[1, 1, :], 1)
 
+        # Index lists not yet implemented
+        with pytest.raises(NotImplementedError, match="list of .* indices"):
+            step_response(
+                sys, timepts=t, input_indices=[0, 1], output_indices=[1])
+
     @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)
     def test_step_response_return(self, tsystem):
         """Verify continuous and discrete time use same return conventions."""
@@ -520,14 +530,24 @@ class TestTimeresp:
         yref = tsystem.yimpulse
         _t, y_00 = impulse_response(sys, T=t, input=0, output=0)
         np.testing.assert_array_almost_equal(y_00, yref, decimal=4)
+
         _t, y_11 = impulse_response(sys, T=t, input=1, output=1)
         np.testing.assert_array_almost_equal(y_11, yref, decimal=4)
 
+        _t, y_01 = impulse_response(
+            sys, T=t, input_indices=[0], output_indices=[1])
+        np.testing.assert_array_almost_equal(y_01, 0 * yref, decimal=4)
+
         yref_notrim = np.zeros((2, len(t)))
         yref_notrim[:1, :] = yref
         _t, yy = impulse_response(sys, T=t, input=0)
         np.testing.assert_array_almost_equal(yy[:,0,:], yref_notrim, decimal=4)
 
+        # Index lists not yet implemented
+        with pytest.raises(NotImplementedError, match="list of .* indices"):
+            impulse_response(
+                sys, timepts=t, input_indices=[0, 1], output_indices=[1])
+
     @pytest.mark.parametrize("tsystem", ["siso_tf1"], indirect=True)
     def test_discrete_time_impulse(self, tsystem):
         # discrete-time impulse sampled version should match cont time
@@ -1396,3 +1416,53 @@ def test_signal_labels():
 
     with pytest.raises(ValueError, match=r"unknown signal name 'x\[2\]'"):
         response.states['x[1]', 'x[2]']         # second index = input name
+
+
+def test_timeresp_aliases():
+    sys = ct.rss(2, 1, 1)
+    timepts = np.linspace(0, 10, 10)
+    resp_long = ct.input_output_response(sys, timepts, 1, initial_state=[1, 1])
+
+    # Positional usage
+    resp_posn = ct.input_output_response(sys, timepts, 1, [1, 1])
+    np.testing.assert_allclose(resp_long.states, resp_posn.states)
+
+    # Aliases
+    resp_short = ct.input_output_response(sys, timepts, 1, X0=[1, 1])
+    np.testing.assert_allclose(resp_long.states, resp_posn.states)
+
+    # Legacy
+    with pytest.warns(PendingDeprecationWarning, match="legacy"):
+        resp_legacy = ct.input_output_response(sys, timepts, 1, x0=[1, 1])
+    np.testing.assert_allclose(resp_long.states, resp_posn.states)
+
+    # Check for multiple values: full keyword and alias
+    with pytest.raises(TypeError, match="multiple"):
+        resp_multiple = ct.input_output_response(
+            sys, timepts, 1, initial_state=[1, 2], X0=[1, 1])
+
+    # Check for multiple values: positional and keyword
+    with pytest.raises(TypeError, match="multiple"):
+        resp_multiple = ct.input_output_response(
+            sys, timepts, 1, [1, 2], initial_state=[1, 1])
+
+    # Check for multiple values: positional and alias
+    with pytest.raises(TypeError, match="multiple"):
+        resp_multiple = ct.input_output_response(
+            sys, timepts, 1, [1, 2], X0=[1, 1])
+
+    # Make sure that LTI functions check for keywords
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        resp = ct.forced_response(sys, timepts, 1, unknown=True)
+
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        resp = ct.impulse_response(sys, timepts, unknown=True)
+
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        resp = ct.initial_response(sys, timepts, [1, 2], unknown=True)
+
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        resp = ct.step_response(sys, timepts, unknown=True)
+
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        info = ct.step_info(sys, timepts, unknown=True)