绘制具有多项式矩阵的系统的响应答案

【问题标题】：Plotting response of a system with polynomial matrices绘制具有多项式矩阵的系统的响应
【发布时间】：2014-03-09 23:23:02
【问题描述】：

当系统、输入、输出矩阵是多项式时，如何绘制系统的时间响应？例如

A(x) = [0.59*x 1.67*x; 0.1 0.2]
B(x) = [2.3*x; 0.3]
C    = [1 0]
Operating region of x = [-2, 2]
Initial condition x(0) = [2 0]

如果矩阵是常数，我可以使用ss 和lsim 来绘制它。但是，在这种情况下我该怎么做呢？我是 Matlab 和控制系统的新手。

【问题讨论】：

您当然可以使用通用 ODE 求解器，例如 ode45 来集成它。或者你可以尝试线性化它。你也试过了吗？
请添加完整的微分方程。似乎您对因变量和自变量都使用了字母 x。是y'(x)=A(x)·y+B(x)，即仍然是线性的，还是y'(x)=A(y)·y+B(y)，可以是任何东西？

标签： matlab control-theory

【解决方案1】：

使用您的示例，我将向您介绍 MatLab 中通用 ODE 模拟的基本思想（这在各种编程语言中非常相似）。

（首先，我假设您在 A 和 B 矩阵中写入的 x 实际上是 x1，因为您没有指定。我还从上下文中假设您的状态向量是 [x1 x2]）。

创建一个函数

当前时间，t_now
当前状态向量，x_now
完整的输入数组，你
全时数组，t

并使用它们来计算状态导数（xdot）。如您所见，仅在索引控制输入 (u) 时才需要全时数组 (t)。基本上，它将是一个编码函数

xdot = f(t, x, u)

这是 ODE 的一般形式。

function xdot = myODE(t_now, x_now, u, t)

    A = [ 0.59*x_now(1), 1.67*x_now(1);
              0.1,            0.2      ] ;

    B = [ 2.3*x_now(1);
             0.3       ] ;

    u_now = interp1(t, u, t_now) ; % get u(t=t_now)

    xdot = A*x_now + B*u_now ;

end

接下来，使用像 MatLab 的 ode45 这样的 ODE 求解器创建一个运行模拟的脚本。如果您想了解这些求解器的工作原理，请阅读数值积分。 MatLab 的 ode45 使用 Runge-Kutta 方法。

%// time range over which to compute simulation
t = [0 : 0.01 : 5] ;

%// define the input signal
%// why not make it zero to look at the natural response
u = zeros(1, length(t)) ;

%// initial condition
x0 = [2, -5] ;

%// call ode45
%// first argument: anonymous call to myODE that tells it which inputs to use
%// second argument: the time range for the simulation
%// third argument: the initial condition
%// first output: the sim time, may be more detailed than the original t
%// second output: the full output including all states
[t_out, response] = ode45(@(t_now, x_now) myODE(t_now, x_now, u, t), t, x0) ;

%// The response contains two column vectors: one for x1 and one for x2 at
%// every time in t. You can extract "the true output" which based on
%// your C matrix is x1.
y = response(:,1) ;

%// Plot results
plot(t, u, 'r--', t_out, y, 'k') ;
grid on ;
legend('input', 'output')
xlabel('time')
ylabel('value')

假设您没有将 u 作为预定义的输入信号，而是将其作为当前状态的函数。然后在 myODE 函数中简单地修改 u_now 的计算。在最简单的情况下，u 根本不是时间的函数，在这种情况下，您甚至根本不需要将 u 或全时数组传递给 myODE！例如，如果

u := -k*x1

那么你的 ODE 函数可以是

function xdot = myODE(t_now, x_now)

    A = [ 0.59*x_now(1), 1.67*x_now(1);
              0.1,            0.2      ] ;

    B = [ 2.3*x_now(1);
             0.3       ] ;

    k = 5 ;
    u_now = -k*x_now(1) ;

    xdot = A*x_now + B*u_now ;

end

调用 ode45 是

[t_out, response] = ode45(myODE(t_now, x_now), t, x0) ;

并且不需要在脚本中定义 u。

希望有帮助！

【讨论】：