[自动控制原理] Matlab求解二阶系统阶跃响应（过阻尼）响应时间ts的ts/T1-T1/T2曲线方法（模拟二阶系统的阶跃响应及其性能分析）

最新推荐文章于 2024-09-23 14:39:52 发布

程序菜鸟一只

最新推荐文章于 2024-09-23 14:39:52 发布

阅读量1.3w

点赞数 7

文章标签： matlab

本文链接：https://blog.csdn.net/sbsbsb666666/article/details/120520274

版权

学习“二阶系统的阶跃响应的性能分析”之后，每次用T1/T2求解ts时，都'只能查表获得'（这是老师说的哈），为了更方便地求得更加精确地ts/T1值，编写了此函数，也对二阶系统的阶跃响应有了一些了解（我可能自控原理不太好QAQ）

最初我考虑直接使用solve来求解，直接输入函数，但matlab提示错误“Unable to find explicit solution”,错误代码及错误界面如下，此时y解不出来

由自控原理的书可以找到ts/T1和T1/T2满足如下关系：

$h(t_{s})=1+\frac{1}{\frac{T_{2}}{T_{1}}-1}e^{-\frac{t_{s}}{T_{1}}}+\frac{1}{\frac{T_{1}}{T_{2}}-1}e^{-\frac{T_1}{T_2}\frac{t_{s}}{T_{1}} }=0.95$ (其中取5%误差)

在“过阻尼”中，T1,T2和阻尼系数 $\xi$ 有如下关系：
$T_1=\frac{1}{\omega_n}\frac{1}{\xi+\surd{(\xi^2-1)}}$

$T2=\frac{1}{\omega_n}\frac{1}{\xi^2+\surd{(\xi^2-1)}}$

故有

$\xi_n=\frac{1}{2\omega_n}(\frac{1}{T_1}+\frac{1}{T_2})$

注意！一定不要把欠阻尼和过阻尼的公式搞混，本文针对过阻尼情形，而对于欠阻尼情形，有公式 $t_s=-\frac{ln0.05+\frac{1}{2}ln(1-\xi ^2) }{\xi \omega_n}\approx \frac{3.5}{\xi \omega_n}$

下面是作者尝试将“过阻尼”的方法代入“欠阻尼”公式，并用 $\xi_n$ 表示出 $\frac{T1}{T2}$ 得到的错误结果(这差不多荒废了作者的一个下午)

后来，也尝试过事先替换值的方法，也失败过很多遍

要注意一点（******matlab中的solve函数不是什么都能解出来的！！！***********）

最终，通过简化方程的方法，将正确的函数曲线画了出来

下面附上我修改了好多遍的函数代码（将这个函数保存为FindResponseTime.m，无输入和输出时，绘制响应曲线，而有输入(输入T1/T2)时，绘制带stem的响应曲线，当有输出(如a=FindResponseTime(3.5))时），计算值并返回

function [result]=FindResponseTime(n)
%输入T1/T2  输出ts/T1
    if nargout==0  %无输出值后，绘制曲线
        syms t
        disp("正在解算并绘制响应函数图像，这可能需要2~3分钟时间.........")
        warning off
        
        x=1.05:0.05:8;
        kexi=(1+x)./(2*sqrt(x));
        coe1=1./(1./x-1);
        coe2=1./(x-1);
        y1=[];
        
        for i=1:1:140
            %----------绘图所需的值矩阵----------
            Eq=coe1(i)*exp(-t)+coe2(i)*exp(x(i))^(-t)==-0.05;
            y1=[y1,solve(Eq,t)];
        end
        
        if nargin==0
        %-------计算图像代码：-----------------------
            hold on
            plot(x,y1,'-b','LineWidth',1.3)
            
            for i=1:1:7
                plot(1+i,y1(20*i),'.r',"MarkerSize",15);
                text(1+i,y1(20*i)+0.2,["\zeta=",num2str(kexi(20*i))]);
            end
            xlabel("T1/T2")
            ylabel("t_s/T_1")
            grid;
            title('二阶系统的过阻尼响应曲线')
            axis([1,8,3,4.8]);
       %--------------------------------
        end
        
        if nargin==1
            disp("请注意:绘制线图时，输入的T2/T1要在1~8之间")
            %-------利用回调计算值-----
            y2=FindResponseTime(n);
            %------------------
            hold on
            plot(x,y1,'-b','LineWidth',1.3)
            for i=1:1:7
                plot(1+i,y1(20*i),'.r',"MarkerSize",15);
                text(1+i,y1(20*i)+0.2,["\zeta=",num2str(kexi(20*i))]);
            end
            stem(n,y2)
            t=(n-1)/0.05
%这里我本想用text和num2str将数据写在旁边，但无法写入，错误代码附在注释中
%text(n,y2-0.05,['t_s/T_1=',num2str(y1(t))])
            
            xlabel("T1/T2")
            ylabel("t_s/T_1")
            grid;
            title('二阶系统的过阻尼响应曲线')
            axis([1,8,3,4.8]);
        end
    end
    
    if nargout==1
    %------------------------
            syms t
            x2=n;
            coff1=1/(1/x2-1);
            coff2=1/(x2-1);
            Eq2=coff1*exp(-t)+coff2*exp(x2)^(-t)==-0.05;
            result=solve(Eq2,t);
    %---------计算输入所相对应的值(输入参数错误时自动提示参数不足)---------
    end
end

这里解决问题的方案是：先单独用x计算出两个e^...前面的系数coe1和coe2,再增加t参数，这样先计算出系数之后，solve的运算量就少了很多（）

后面的coff是避免重名的

最后说一下使用中的注意事项

对于FindResponseTime(k) 要求k应当在1-8之间（不然会超出图去）

但是a=FindResponseTime(k),你可以随便输入大于1的（如a=FindResponseTime(9)也没问题）（另外这个比其他两个命令运行快）

加上最终成功的运行截图：