1、系统的思想
系统的思想从古到今一直发生着演变。中国古代文化中的系统思想主要体现在对世界的认识和解释上，如《道德经》中的“道生一，一生二，二生三，三生万物”，论证了世间万物由要素组成，要素之间存在相互联系相互作用；又如春秋战国诸子百家对军事、自然、社会等现象进行讨论，形成了中国独有的军事、自然科学和社会伦理体系，并将其应用于具体的工程项目。近现代西方的系统思想则更关注系统的运行规律，其数学、物理、化学、生物学等基础科学和自然科学的发展，使人类对自然系统的组成、要素关系和运行规律等方面的认识取得了突破性成果，发展和完善了科学的方法论。
2、系统的定义
这里引用钱学森的观点给出“系统”的定义：系统是由相互作用相互依赖的若干部分结合而成的、具有特定功能的有机整体，而且这个有机整体又是它从属的更大系统的组成部分。由此可以提炼出认识一个系统的关键点在于组成系统的要素和各要素如何相互作用。根据系统层次结构的复杂程度，可以将系统分为简单系统和复杂系统：直接由要素构成的两层次系统称为简单系统，由多层子系统按照一定结构构成的系统称为复杂系统。例如，中国载人航天工程就是一个复杂系统，由航天员系统、空间应用系统、载人飞船系统、货运飞船系统、长征二号F运载火箭系统、长征七号运载火箭系统、长征五号B运载火箭系统、酒泉发射场系统、文昌发射场系统、测控通信系统、空间实验室系统、空间站系统、着陆场系统和光学舱系统等十四大系统组成，各个系统又由若干个子系统构成，如测控通信系统就包含雷达系统、光学系统等子系统。
图1 中国载人航天工程系统组成
此外，按照时间和状态变量的取值特性，可以将系统分为连续、离散系统和混合系统状态变量随时间连续变化的系统称为连续系统；时间变量和状态变量均只取有限个离散值的系统称为离散系统；兼具连续和离散成分的系统称为混合系统。按照系统要素之间的关系或模型的形式，可以将系统分为线性系统和非线性系统：要素之间的关系能够用线性数学模型描述的系统称为线性系统，否则称为非线性系统。除此之外，还有其他的对系统的分类方式，如封闭系统和开放系统、动态系统和静态系统等。
[1]中国载人航天官方网站 中国载人航天工程网 http://www.cmse.gov.cn/gygc/xtzc/htyxt/。
[2] 钱学森．论宏观建筑与与微观建筑：杭州出版社，2001
为了研究、分析、设计和实现一个系统，必然需要进行对系统的试验。试验的方法可以分为两大类：一种是直接在真实系统上进行试验，一种则是根据系统的组成要素和要素间的相互作用，将系统抽象为模型，对模型进行试验。系统模型通常分为物理模型、概念模型和数学模型。本文重点论述基于数学模型描述的系统。
通过对系统建模进行试验和结果分析，就是所谓的“系统仿真”。直接在真实系统上进行试验，显然更有助于系统的研究分析，而系统仿真不可能完全符合系统实际运行情况，那么为什么要进行系统仿真呢？一言以蔽之，系统仿真为系统设计和系统分析提供了成本较低的有力支持。
以载人航天工程为例，其一，在系统设计阶段，真实的系统尚未建立，如何了解和验证系统的性能？载人航天作为举国体制下的重要工程，牵一发而动全身，在给出最终系统方案之前，必然要先确定所采用方案中系统的性能，但若为每一个方案都搭建一套真实的系统，在资源和人力上显然是不可承受的，系统仿真在这种情况下就为系统设计提供了有力支持；其二，在真实系统上进行试验可能引起系统破坏或故障，如果用真实的飞船搭载火箭进行发射试验，过程中不可避免会对系统产生破坏，其成本必然是高昂的，通过系统仿真可以减少试验成本；其三，在实际的发射过程会遇到各种复杂的随机问题，如何在试验阶段对其进行预测和提出解决方案呢？另一方面，某些系统在多次实验中要求相同的试验条件，如此才能更好地评估性能和发现问题，这在实际环境下是不可能实现的，而系统仿真则可以对同一试验条件的精准再现；其四，对于载人航天来说，每一次发射都要耗费漫长的时间和高昂的金钱成本，次次都在真实的系统上进行试验成本过高，系统仿真既节省了试验时间，又降低了试验成本。
仿真的基本方法是建立系统的结构模型和量化分析模型，对模型进行仿真试验。根据模型的不同，可以将系统仿真方法分为三类：物理仿真、数学仿真和半实物仿真。物理仿真是按照真实系统的物理性质构造系统的物理模型，并在物理模型上进行试验。物理仿真虽然直观形象，但具有较多实验限制，且成本较高。数学仿真是对实际系统进行抽象，并将其特性用数学关系加以描述而得到系统的数学模型，将其转换为适合计算机编程的仿真模型，进行仿真试验。现代计算机技术的发展为系统仿真提供了条件。数学仿真在通信、控制、电力等多个工程和科研领域具有广泛应用。图1、图2是系统仿真的简单示例。
图1 双质量-弹簧-阻尼系统
图2 低通RC滤波器
由于连续系统和离散系统的数学模型具有很大差别，所以数学仿真可以分为连续系统仿真和离散系统仿真。半实物仿真则是对物理仿真和数学仿真的综合应用，即对系统中运行规律明确的部分且简单的部分进行数学仿真，规律不明确或难以数学建模的部分进行物理建模，二者联合起来完成仿真试验。例如在通信领域，对一个通信设备的系统仿真，除了通过数学建模完成对其通信系统的仿真分析，还要通过对设备使用环境进行物理仿真，以确保设备在真实环境下的运转情况。
系统仿真是科研和工程领域不可缺少的环节，对工程技术和科学研究的发展具有重要意义。
[1] 福州大学《系统仿真技术》研究生课程
[2] http://www.bilibili.com/video/BV16E411A7zv?p=3&vd_source=f52c61024bca044c638d9c2fb92b3ede.
系统仿真软件是一种用于模拟实际系统的软件，它可以帮助我们更好地理解系统的运行情况。仿真软件是用于建模和分析复杂系统的强大工具。它被用于从工程到金融的各种行业，以帮助了解系统如何工作以及如何改进。不同的仿真软件有不同的特点，因此在选择仿真软件时，我们需要考虑不同的系统仿真软件的优缺点。
MATLAB/Simulink是一款功能强大的系统仿真软件，它可以帮助我们快速构建复杂的系统模型，并且可以模拟多种不同的系统。MATLAB是一个功能强大的通用工具，可用于模拟各种系统。它易于使用，具有广泛的功能，包括图形用户界面、强大的编程语言和预构建模型库。MATLAB也是高度可扩展的，允许用户创建自己的模型和仿真。此外，MATLAB/Simulink还提供了丰富的可视化工具，可以帮助我们更好地理解系统的运行情况。Simulink是一个图形编程环境，允许用户创建和模拟复杂系统。它被设计为易于使用，并具有广泛的功能，包括一个预先构建的模型库、一种强大的编程语言和一个图形用户界面。Simulink还具有高度的可扩展性，允许用户创建自己的模型和仿真。但是，MATLAB/Simulink的上手较难，对于初学者来说，可能需要花费更多的时间来学习。
LabVIEW是另一款功能强大的系统仿真软件，它可以帮助我们快速构建复杂的系统模型，并且可以模拟多种不同的系统。此外，LabVIEW还提供了丰富的可视化工具，可以帮助我们更好地理解系统的运行情况。LabVIEW主要用于开发测量或控制系统。
Modelica是一种用于建模和仿真复杂系统的面向对象的编程语言。它提供了一种灵活的方法来描述系统的行为，并且可以用于模拟和分析复杂的系统。Modelica语言支持多种类型的模型，包括动态系统、静态系统、热系统、流体系统和电气系统。它还支持多种类型的仿真，包括时间域仿真、频率域仿真和状态空间仿真。
OpenModelica属于开源模拟软件包。OpenModelica是一个功能强大的多功能工具，可用于模拟各种系统。它被设计为易于使用，并具有广泛的功能，包括一个预先构建的模型库、一种强大的编程语言和一个图形用户界面。
Dymola基于Modelica语言提供一个基于模型的仿真环境，它可以用于模拟和分析复杂的系统。它提供了一个可视化的编程环境，可以帮助用户快速构建模型，并使用图形化的工具进行模拟和分析。与Matlab/Simulink不同，Dymola支持模型驱动的仿真，可以更好地模拟复杂的系统。此外，Dymola还支持多种仿真技术，如动态系统建模、模型验证和可视化等。
Scilab是一款开源的科学计算软件，它可以用于数值计算、矩阵运算、绘图、科学可视化、数据分析、模拟仿真等。它拥有丰富的函数库，可以满足各种科学计算需求。Scilab可以运行在Windows、Linux、Mac OS X等操作系统上，支持多种编程语言，如C、C++、Fortran等，可以与其他软件进行交互，如Matlab、Maple、Mathematica等。
总体来说，对于控制、通信、电气等相关领域的理论研究推荐MATLAB/Simulink，LabVIEW则推荐用于测控系统，Modelica是一种编程语言，其自带的求解器较为弱，偏好该语言的用户可以尝试Dymola等相关软件，Scilab则是一个较为小众的开源软件。
一、 最小二乘原理介绍
最小二乘法是一种常见的数学方法，用于解决线性回归问题。它可以用来估计因变量（目标变量）与自变量之间的线性关系，即通过已知的一组数据来确定未知的回归方程。
具体来说，最小二乘法的目标是找到一条直线（或者更一般的曲线），使得这条直线与已知数据点的距离平方和最小。这条直线的方程称为回归方程，可以用来预测因变量的值。
最小二乘法的原理是通过最小化误差平方和来确定回归方程的系数。误差指的是已知数据点到回归直线的距离，也就是因变量的真实值与预测值之间的差异。通过最小化误差平方和，可以找到最优的回归系数，使得预测值与真实值之间的误差最小。
在实际应用中，最小二乘法常常用于数据拟合和预测。例如，在金融领域，最小二乘法可以用来预测股票价格和交易量；在工程领域，最小二乘法可以用来拟合实验数据和优化设计参数。
在计算中，最小二乘法可以通过矩阵运算来求解。具体来说，可以使用线性代数中的矩阵求解方法，例如求解矩阵的逆或者使用QR分解。在北太天元中，可以使用polyfit函数来实现最小二乘法求解线性回归问题。
二、 根据最小二乘原理进行多项式拟合
原始数据点为，根据最小二乘法进行多项式拟合：
(1) 设拟合多项式为
(2) 各点到这条曲线的距离之和，即误差平方和为：
(3) 根据最小二乘原理，等式右边求ai偏导数：
(4) 化简得：
(5) 化成矩阵形式：
(6) 化简得
(7) 求解上述方程组，即可得到系数矩阵，从而得到拟合多项式。
三、 基于北太天元实现最小二乘法，并对比拟合结果与polyfit结果
附录代码：
clear all
close all
clc
N=10;%多项式阶数
%原始数据点
X=10:0.05:15;
M=length(X);
p=[5.6,2.5,3.3,0.18,0.29,2.4,4.8,1.2,7.1,1];
Y=zeros(1,M);
for i=1:10
Y=Y+p(i)*X.^(i-1);
end
% Y=Y+unifrnd(-2e9,2e9,1,M);
%多项式拟合
X1=zeros(N+1,M);
for i=1:M
X1(1,i)=1;
end
for i=2:N+1
X1(i,:)=power(X,i-1); %构造系数矩阵
end
XX=X1';
P_n=(XX'*XX)\XX'*Y';
X_n=min(X):max(X);
Y_n=zeros(1,length(X_n));
for i=1:N+1
Y_n=Y_n+P_n(i)*X_n.^(i-1);
end
%polyfit
P=zeros(1,N);
P=polyfit(X,Y,N);
X_polyfit=zeros(1,M);
Y_polyfit=zeros(1,M);
X_polyfit=linspace(min(X),max(X),M);
Y_polyfit=polyval(P,X_polyfit);
figure(1);
plot(X,Y,'.');
hold on
plot(X_n,Y_n,'g');
hold on
plot(X_polyfit,Y_polyfit,'r');
legend('原始数据点','多项式拟合','polyfit','Location','northwest');
该求解器有变步长和定步长两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法。
ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长（变步长）的常微分方程数值解法，其整体截断误差为(Δx)^5。解决的是非刚性常微分方程。
ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode15s试试。
在各种龙格－库塔法当中有一个方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。
令初值问题表述如下。
y' =f(t,y)，y(t0)=y0
则，对于该问题的RK4由如下方程给出：
y(n+1)=y(n)+h/6*(k1+2*k2+2*k3+k4)
其中
k1=f(t(n),y(n))
k2=f(t(n)+h/2,y(n)+h/2*k1)
k3=f(t(n)+h/2,y(n)+h/2*k2)
k4=f(t(n)+h,y(n)+h*k3)
这样，下一个值（yn+1）由现在的值（yn）加上时间间隔（h）和一个估算的斜率的乘积所决定。该斜率是以下斜率的加权平均：
k1是时间段开始时的斜率；
k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；
k3也是中点的斜率，但是这次采用斜率k2决定y值；
k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值：
RK4法是四阶方法，也就是说每步的误差是h阶，而总积累误差为h阶。
求解微分方程d²y/dt²- 7(1-y²)dy/dt + y = 0在初始条件下y(0) = 1、y’(0) = 0下的解，并画出解的图。
解：设x1 = y，x2 = dy/dt，将二阶方程化成一阶方程组
dx1 = x2 x1(0) = 1;
dx1/dt = 7(1-x1²)*x2-x1 x2(0) = 0;
在北太天元脚本编辑器窗口输入以下程序，并保存文件名为vdp.m。
function fy = vdp(t,x)
fy = [x(2);7*(1-x(1)^2)*x(2)-x(1)];
在北太天元命令行窗口输入以下程序。
>> y0 = [1;0];
>> [t,x] = ode45(@vdp,[0 40],y0);
>> y = x(:,1);
>> plot(t,y)
>> xlabel('t')
>> ylabel('y')
输出结果如下图所示：
用雷达跟踪目标，目标的运动可以看成是在径向和横向内的二维运动，其运动方程和观测方程分别为：
s1(k)、v1(k)和y1(k)分别为径向距离、速度和观测值，而s2(k)、v2(k)和y2(k)分别为横向距离、速度和观测值。u1(k)和u2(k)是状态噪声，是目标速度的波动；w1(k)和w2(k)是观测噪声；四种噪声的均值都为0，呈高斯分布，互不相关。
状态方程和输出方程：
其中：
在推导上式时，wi(k)和wi(k+1)是随机过程中不同时刻的两个随机变量，我们认为这两个随机变量统计独立，而且w(k)是平稳随机过程，其不同时刻的方差相同。两个时刻的时差T为雷达扫描一圈的时间。
仿真结果如下：
clear allclose allclcT=1;sigma_u1=300;sigma_u2=0.12;sigma_w1=900000;sigma_w2=900000;N=500;u1=sqrt(sigma_u1).*randn(1,N);u2=sqrt(sigma_u2).*randn(1,N);w1=sqrt(sigma_w1).*randn(1,N);w2=sqrt(sigma_w2).*randn(1,N);A=[ 1,T,0,0; 0,1,0,0; 0,0,1,T; 0,0,0,1; ];C=[ 1,0,0,0; 0,0,1,0; ];Q=[ %状态噪声协方差矩阵 0,0,0,0; 0,sigma_u1,0,0; 0,0,0,0; 0,0,0,sigma_u2; ];R=[ %观测噪声协方差矩阵 sigma_w1,0; 0,sigma_w2; ];I=[ 1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0,1; ];x0=[10000;300;900;256;];%初始值% y0=[10000;0;];x=zeros(4,N);y=zeros(2,N);x(:,1)=x0;% y(:,1)=y0;x_calculation=zeros(4,N); %状态滤波估计值u=[zeros(1,N);u1;zeros(1,N);u2];w=[w1;w2];y(:,1)=C*x(:,1)+w(:,1);P2=[ sigma_w1, sigma_w1/T, 0, 0; sigma_w1/T,sigma_u1+2*sigma_w1/(T^2),0, 0; 0, 0, sigma_w2, sigma_w2/T; 0, 0, sigma_w2/T,sigma_u2+2*sigma_w2/(T^2); ];p=P2;for i=1:N-1 x(:,i+1)=A*x(:,i)+u(:,i);%真实值 y(:,i+1)=C*x(:,i+1)+w(:,i+1); endx_calculation(:,2)=[y(1,1) , (y(1,2)-y(1,1))/T ,y(2,2), (y(2,2)-y(2,1))/T]';for i=3:N x1=A*x_calculation(:,i-1);%状态预测方程 p1=A*p*A'+Q; %状态预测误差协方差矩阵 K=p1*C'*inv(C*p1*C'+R); %最佳增益方程 x_calculation(:,i)=x1+K*(y(:,i)-C*x1); %滤波估值方程 p=(I-K*C)*p1; %滤波估值误差 协方差方程end t=1:N; figure(4)plot(t,y(1,:),'b',t,x_calculation(1,:),'r');xlabel('t');ylabel('s1');legend('径向距离观测值','径向距离估计值','Location','southeast');title('二维卡尔曼滤波器');figure(5)plot(t,y(2,:),'b',t,x_calculation(3,:),'r');xlabel('t');ylabel('s2');legend('横向距离观测值','横向距离估计值','Location','southeast');1) 忽略转向系的影响，以前、后轮转角作为输入；
2) 汽车只进行平行于地面的平面运动，而忽略悬架的作用；
3) 汽车前进（纵轴）速度不变，只有沿y轴的侧向速度和绕z轴的横摆运动(ay<0.4g) ；
4) 驱动力不大，对侧偏特性无影响；
5) 忽略空气阻力；
6) 忽略左右轮胎因载荷变化引起轮胎特性的变化；
7) 忽略回正力矩的变化。
根据模型假设建立如图1所示的二自由度汽车模型。
对模型受力分析，存在3个方向的受力平衡，分别为x、y和绕Z的力矩平衡，建立力学方程如下。
在baltamulink中搭建状态空间模型，模型如图所示。
（1）在前轮偏转角为1，后轮偏转角为1，车速为40km/h的情况下，输出前后轮的横向位移情况，输出结果如图3。
图3
通过建立汽车动力学模型，对汽车操纵性进行饿模拟。根据仿真结果可以发现车速和前轮转角都对二自由度汽车的操纵稳定性有很大影响。汽车以较低速度、较小的前轮转角行驶时，是相对安全的。
通过分析图3可以看出前、后轮的横向位移都是发散的，这是因为给前轮的一个阶跃响应，一直存在前轮转角，同时系统没有加入闭环控制，属于开环控制，这就导致前后轮的横向位移处于发散状态。
baltamatica代码如下
clc;clear;close all;
%% 基本车辆参数
v=40/3.6;%输入为km/h，方程单位为m/s
m=16000;%车重
I=10.85*m;%转动惯量
cf=340000;%侧偏刚度
cr=cf;
lf=2.65;%前轴到重心的距离
lr=3.35;%后轴到重心的距离
a=pi/180;
%% 组装矩阵
A=[-(cf+cr)/(m*v) -1-(cf*lf-cr*lr)/(m*v^2) 0 0
-(cf*lf-cr*lr)/I -(cf*lf^2+cr*lr^2)/(I*v) 0 0
v lf 0 0
v -lr 0 0];
B=[cf/(m*v) cr/(m*v)
cf*lf/I cr*lr/I
0 0
0 0];
C=[0 0 1 0
0 0 0 1];
D=zeros(2,2);
问题：单质量弹簧阻尼机械系统类似于下图，外力f(x)为输入量；质量位移x(t)为输出量；质量m为1kg；刚度为5N/m；阻尼系数f为0.3N·s/m。绘制系统位移输出响应曲线。
图1：单质量弹簧阻尼机械系统
首先：单质量弹簧阻尼机械系统的震动方程为：
m dx²/dt² + c dx/dt + kx = f(x)
取状态向量为X(t) = [x(t) dx/dt]’，输出向量为U=[f(x)]，输出向量为Y=[x(t)]，则状态方程为：
dX/dt = AX(t) + BU
Y = CX(t) + DU
式中，A = [0 1; -k/m -c/m]; B = [0; 1/m]; C = [1 0]; D = 0;(m = 1; k = 5; c = 0.3)。
其次，通过计算得到传递函数的分子分母系数为：[1]； [1 0.3 5];
即，传递函数为
G(s) = 1/(s^2 + 0.3 s + 5);
最后，利用北太真元建立传递函数仿真模型，如下图所示：
设置：仿真时长：10s
求解器为：定步长 ode4；
步长：1s
输入常量10
仿真结果如下图所示：
在MATLAB中创建相同模型，仿真时间、仿真求解器、步长均相同；仿真模型和仿真结果如下图所示：
通过对比发现，北太真元计算结果与MATLAB仿真结果完全一致。
问题：
二阶电路动态系统中，该系统是由电阻R、电感L和电容C组成的无源网络；ui(t)为输入，i(t)为电流，u0(t)为输出；设R=1Ω，L= 2H，C=2F；系统的初始状态为0，外加的输入为单位阶跃信号。求系统的输出波形。
图：二阶电路动态系统
首先，在解决该问题时，需要了解基尔霍夫电压定律；它的内容是，在任何一个闭合回路中，各元件上的电压降的代数和等于电动势的代数和，即从一点出发绕回路一周回到该点时，各段电压的代数和恒等于零，即∑U=0。
基尔霍夫电压定律表明：
沿着闭合回路所有元件两端的电势差（电压）的代数和等于零。
或者描述为：
沿着闭合回路的所有电动势的代数和等于所有电压降的代数和。
以方程表达，对于电路的任意闭合回路有：
其中，m 是这闭合回路的元件数目，vk 是元件两端的电压，可以是实数或复数。
基尔霍夫电压定律不仅应用于闭合回路，也可以把它推广应用于回路的部分电路。
根据上述原理，可写出上面问题的回路方程如下：
L di(t)/dt + 1/C∫i(t)dt + Ri(t) = ui(t)
消去中间变量i(t)，可以得到描述网络输入输出关系的微分方程为：
LC d²u0(t)/dt² + RC du0(t)/dt + u0(t) = ui(t)
带入电路参数R=1Ω，L= 2H，C=2F；整理可得：
d²u0(t)/dt² + 0.5 du0(t)/dt + 0.25u0(t) = 0.25ui(t)
从方面方程可以得到，传递函数分子等式右边系数：[0.25]；分母为左边等式系数：[1 0.5 0.25]
根据该系数，在北太真元设计传递函数仿真模型，再加入一个阶跃信号模块，得到仿真模型如下图所示：
设置：仿真时长：40s
求解器为：定步长；
步长：1s；
得到仿真结果如下图所示：
在MATLAB中创建相同模型，仿真时间、仿真求解器、步长均相同；仿真模型和仿真结果如下图所示：
通过对比发现，北太真元计算结果与MATLAB仿真结果完全一致。
由于汽车在行走时，路面不平，汽车行驶中的路面可简化成正弦函数。
可把汽车行走的路面看做激励，忽略轮胎的弹性与质量，得到分析车身垂直振动的最简单的单质量系统，适用于低频激励情况。汽车行驶可看作如下模型：
上图为单一自由度系统的简图，设X(t)及Xs(t)分别是质量块及支承的位移，支承的运动规律是：
Xs =asinwt
由于支承的运动，质量块收到的弹性恢复力为k(X - Xs)，阻尼力为c(Vx-Vxs)
根据达朗伯原理可得如下的运动微分方程：
由（1）和（2）得：
在此系统中除了有弹性恢复力及阻尼力作用外，还始终作用于简谐激励力：
Px=P0sinwt
简谐激励：
激励随时间的变化规律可用正弦或余弦函数表示；
振动响应亦为时间的正弦和余弦函数（简谐振动）。
结合上面的运动微分方程和简谐激励力方程，可得系统的运动微分方程为：
令：物体质：m = 1 kg，弹簧刚度：k = 3 N/m，阻尼：c = 4 N·s/m，作用力P = 2sin(2t + π/3)，研究物体的位移随时间的变化规律。
通过北太真元建立系统运动微分方程模型，如下图所示：
设置参数：
正弦波产生模块：幅值：2；偏置：0；频率(弧度/秒)：2；相位(弧度)：pi/3≈1；
一阶积分模块：积分初始值：0.5；
增益模块：增益数值：4；
常量模块：常量值：3；
结束时间：10s；
求解器：ode4；
步长：0.01s。
得到的仿真结果，如下图所示：
问题：利用汽车横摆角速度传递函数和质心侧偏角传递函数，对汽车时域响应进行仿真，绘制汽车横摆角速度和质心偏侧角的时域特性曲线。汽车时域响应仿真所需参数见下表。
由于汽车横摆角速度传递函数（G1(s)）和质心偏侧角传递函数（G2(s)）分别为：
G1(s) = ((s - a11)*b21 + a21*b11) / s²- (a11 + a22)*s + a11*a22 - a12*a21;
G2(s) = ((s - a22)*b11 + a12*b21) / s²- (a11 + a22)*s + a11*a22 - a12*a21;
式中，
a11 = (Ka1 + Ka2) / mu; a12 = (aKa1 - bKa2 - mu²) / mu²;
a21 = (aKa1 - bKa2) / Iz; a22 = (a²Ka1 + b²Ka2) / Iz*u;
b11 = -Ka1 / mu; b21 = -aKa1/Iz；
汽车速度分别选取10m/s、20m/s、30m/s；在仿真时间0s时给前轮一个阶跃信号，使前轮转角从0°转到15°，并保持不变。根据汽车横摆角速度传递函数和质心偏移角传递函数，建立模型，绘制不同车速下的汽车横摆角度和质心侧偏角的时域特性曲线。
首先：通过北太天元计算汽车横摆角度速度传递函数分子和分母的系数，在北太天元依次输入下面语句；
>> m=3018; Iz=10437; a=1.84; b=1.88; k1=-23147; k2=-38318;
>> u= [10 20 30];
>> a11 = (k1 + k2)/m./u;
>> a12 = (a*k1 - b*k2 -m.*u^2)/m./u^2;
>> a21 = (a * k1 - b * k2)/Iz;
>> a22=(a^2*k1 + b^2*k2)/Iz./u;
>> b11 = -k1/m./u;
>> b21 = -a*k1/Iz;
>> b1 = b21;
>> b2 = a21*b11-a11*b21;
>> b3 = -a11-a22;
>> b4 = a11.*a22-a12.*a21;
>> num = [b1 b2];
>> den = [1,b3,b4];
得到结果如下图1所示 ；将命令行窗口，和工作区窗口放大后如图2、图3所示；
图1
图2
图3
因为，汽车横摆角度速度传递函数的分子系数为：num = [b1, b2]; 分母系数为：den = [1, b3, b4]; 所以，从图3红色框中可以得到各项系数如下：
当汽车速度 s = 10 m/s 时，分子系数：num = [4.0807 10.4748]; 分母系数：den = [1 4.0851 6.7181];
当汽车速度 s = 20 m/s 时， 分子系数：num = [4.0807 5.2347]; 分母系数：den = [1 2.0425 3.7956];
当汽车速度 s = 30 m/s 时， 分子系数：num = [4.0807 3.4916]; 分母系数：den = [1 1.3617 3.2544];
又因为，在仿真时间0s时给前轮一个阶跃信号，使前轮转角从0°转到15°；所以模型还需一个阶跃信号模块，阶跃时间=0；且，还需一个增益模块，增益= pi*15/180 = 0.2618。
通过北太真元建立汽车横摆角度速度传递函数模型，如下图所示：
设置参数：
仿真时长：10s；步长0.01s；求解器：ode4
得到的仿真结果，如下图所示：
紫色代表速度10m/s时响应曲线；
墨绿色代表速度20m/s时响应曲线；
橙色代表速度30m/s时响应曲线；
首先：通过北太天元计算质心偏侧角传递函数分子和分母的系数，在北太天元依次输入下面语句；
>> m=3018;Iz=10437;a=1.84;b=1.88;k1=-23147;k2=-38318;
>> u= [10 20 30];
>> a11 = (k1 + k2)/m./u;
>> a12 = (a*k1 - b*k2 -m.*u.^2)/m./u.^2;
>> a21 = (a * k1 - b * k2)/Iz;
>> a22=(a^2*k1 + b^2*k2)/Iz./u;
>> b11 = -k1/m./u;
>> b21 = -a*k1/Iz;
>> b1 = b11;
>> b2 = a12*b21-a22*b11;
>> b3 = -a11-a22;
>> b4 = a11.*a22-a12.*a21;
>> num = [b1 b2];
>> den = [1,b3,b4];
命令行窗口如图4所示；参数计算结果如图5所示；
图4
图5
因为，质心偏侧角传递函数的分子系数为：num = [b1, b2]; 分母系数为：den = [1, b3, b4]; 所以，从图5红色框中可以得到各项系数如下：
当汽车速度 s = 10 m/s 时，分子系数：num = [0.766965 -2.11146]; 分母系数：den = [1 4.08507 6.71806];
当汽车速度 s = 20 m/s 时， 分子系数：num = [0.383482 -3.58841]; 分母系数：den = [1 2.04254 3.7956];
当汽车速度 s = 30 m/s 时， 分子系数：num = [0.255655 -3.86191]; 分母系数：den = [1 1.36169 3.2544];
又因为，在仿真时间0s时给前轮一个阶跃信号，使前轮转角从0°转到15°；所以模型还需一个阶跃信号模块，阶跃时间=0；且，还需一个增益模块，增益= pi*15/180 = 0.2618。
通过北太真元建立质心偏侧角传递函数模型 同 汽车横摆角度速度传递函数模型。
设置参数：
仿真时长：10s；步长0.01s；求解器：ode4
得到的仿真结果，如下图所示：
浅绿色表速度10m/s时响应曲线；
紫代表速度20m/s时响应曲线；
墨绿色代表速度30m/s时响应曲线；
实例：已知传递函数G(s) = ω²/ s²+ 2ζωs + ω²，分析阻尼系数和固有频率对性能的影响。
（1）假设 ω = 1, ζ = 0, 0.8, 1.5;
（2）假设 ζ = 1, ω = 1, 2, 3;
从（1）可得，阻尼系数传递函数的系数可以是：
G(s1) = [1; 1 0 1];
G(s2) = [1; 1 1.6 1];
G(s3) = [1; 1 3 1];
通过北太真元建立“阻尼系数对系统性能的影响”模型，如下图所示：
设置参数：仿真时长：20s；步长0.1s；求解器：ode4
得到的仿真结果，如下图所示：
墨绿色代表阻尼系数=0时的特性曲线；
淡绿色代表阻尼系数=0.8时的特性曲线；
红色代表阻尼系数=1.5时的特性曲线。
从结果图中可以看出，阻尼系数决定了系统的振荡幅度，阻尼系数越小，振荡幅度越大。
从（2）可得，固有频率传递函数的系数可以是：
G(s1) = [1; 1 0.5 1];
G(s2) = [4; 1 1 4];
G(s3) = [9; 1 1.5 9];
通过北太真元建立“固有频率对系统性能的影响”模型，同上模型；
设置参数：仿真时长：20s；步长0.1s；求解器：ode4
得到的仿真结果，如下图所示：
墨绿色代表固有频率=1时的特性曲线；
淡绿色代表固有频率=2时的特性曲线；
紫色代表固有频率=3时的特性曲线。
从结果图中可以看出，固有频率决定了系统的振荡频率，固有频率越大，系统的振荡越高，响应速度也越快。
问题：利用汽车横摆角速度传递函数和质心侧偏角传递函数，对汽车时域响应进行仿真，绘制汽车横摆角速度和质心偏侧角的时域特性曲线。汽车时域响应仿真所需参数见下表。
取状态向量为X = [β ωr]’，输入向量U = [δ1]，输出向量为Y = [β ωr]’，状态空间方程为：
式中，A = [(K1+K2)/mu, (aK1 - BK2)/mu²-1; (aK1-bK2)/Iz, (a²K1+b²K2)/Iz*u] 称为系统矩阵；B = [-K1/mu; -aK1/Iz] 称为控制矩阵；C = [1 0; 0 1] 称为输出矩阵；D = [0; 0] 称为传递矩阵。
汽车速度分别选取20m/s、30m/s、40m/s；在仿真时间0s时给前轮一个阶跃信号，使前轮转角从0°转到10°，并保持不变。根据汽车状态空间模型，建立模型，绘制不同车速下的汽车横摆角速度和质心侧偏角的时域特性曲线。
首先：通过北太天元计算汽车状态空间方程的系统矩阵和控制矩阵，在北太天元依次输入下面语句；
>> m=2050;Iz=5600;a=1.5;b=1.8;L=3.3;
>> k1=-38900;k2=-39200;
>> u= [20 30 40];
>> a11 = (k1 + k2)/m./u;
>> a12 = (a*k1 - b*k2 -m.*u.^2)/m./u.^2;
>> a21 = (a*k1 - b*k2)/Iz;
>> a22=(a^2*k1 + b^2*k2)/Iz./u;
>> b11 = -k1/m./u;
>> b21 = -a*k1/Iz;
得到结果如下图1所示 ；
图1
将命令行窗口，和工作区窗口放大后如图2、图3所示；
图2
图3
因为，汽车状态空间方程的系统矩阵为：A = [a11, a12; a21, a22]；控制矩阵为：B = [ b11; b21]； 所以，从图3红色框中可以得到各项系数如下：
当汽车速度 s = 20 m/s 时，系统矩阵：A = [-1.90488,-0.98511;2.18036,-1.91547]; 控制矩阵：B = [0.94878; 10.4196];
当汽车速度 s = 30 m/s 时，系统矩阵：A = [-1.26992,-0.993382;2.180436,-1.27698]; 控制矩阵：B = [0.63252; 10.4196];
当汽车速度 s = 40 m/s 时，系统矩阵：A = [-0.952439,-0.996277;2.180436,-0.957737]; 控制矩阵：B = [0.47439; 10.4196];
状态方程输出矩阵C = [1 0; 0 1]；传递矩阵D = [0; 0]。
又因为，在仿真时间0s时给前轮一个阶跃信号，使前轮转角从0°转到10°；所以模型还需一个阶跃信号模块，阶跃时间=0；且，还需一个增益模块，增益= pi*10/180 = 0.1745。
通过北太真元建立汽车状态空间模型，如下图所示：
设置参数：仿真时长：10s；步长0.01s；求解器：ode4
得到的仿真结果，如下图所示：
上半部分代表汽车横摆角速度时域特性曲线；即：
墨绿色代表速度20m/s时的特性曲线；
绿色代表速度30m/s时的特性曲线；
红色代表速度40m/s时的特性曲线。
下半部分代表汽车质心侧偏角时域特性曲线；即：
紫色代表速度20m/s时的特性曲线；
橙色代表速度30m/s时的特性曲线；
蓝色代表速度40m/s时的特性曲线。
使用4:1衰减曲线法设计下列被控传递函数的PI控制器，分别计算P控制、PI控制的参数值，并绘制控制前后系统的单位阶跃响应曲线。4:1衰减法控制参数计算公式如下表所示：
4:1衰减法控制
被控传递函数方程如下：
Gp(s) =1 / 100^3 + 80^s + 17s + -1;
调节参数时，比例系数由小变大，并增加扰动观察响应过程，知道响应曲线峰值衰减比为4:1，记录此时的比例系数Kp为Ks，两个峰值之间的时间周期为周期Ts。
假设：响应曲线峰值衰减比为4:1时的比例系数 Ks= 4.74;
两个峰值之间的时间周期 Ts = 21.9967；
则，按照上面4:1衰减法控制表中的计算公式可得：
P控制：比例系数 Kp = Ks = 4.74；
PI控制：比例系数 Kp = Ks/1.2 = 3.95；
积分时间常数 Ti = 0.5 * Ts = 10.9984;
积分系数 Ki = Kp / Ti = 3.95 / 10.9984 = 0.3591;
P（比例）控制
首先，把控制器设置成纯比例控制，即令积分系数Ki和微分系数Kd为零，在北太真元建立模型，形成比例控制系统，结构如下图所示；
设置参数：
仿真时长：30s；步长0.01s；求解器：ode4
得到的仿真结果，如下图所示：
PI（比例积分）控制
首先，在纯比例控制系统的基础上增加积分系数Ki，令微分系数Kd为零，在北太真元建立模型，形成比例控制系统，结构如下图所示；
设置参数：
仿真时长：30s；步长0.01s；求解器：ode4
得到的仿真结果，如下图所示：
公众号