系统仿真概述

系统与模型

系统

系统是指物质世界中既相互制约又相互联系着的,能够实现某种目的的一个整体

也就是说系统是一个由多个部分组成的,按一定规律连接的,具有特定功能的整体

任何系统都存在四个方面的内容:实体,属性,活动和环境

模型

系统模型实质上是一个由研究目的所确定的,关于系统某一方面本质属性的抽象和简化,并以某种表达形式来描述

系统模型一般可以分为:物理模型和数学模型

为了分析研究一个系统,通常会进行各种实验,其一是在已经存在的真实系统上进行,其二是通过构造模型,利用模型试验的方法进行

其中,第二种形式的实验越来越多,理由如下:

  1. 系统还处于设计阶段时,真实系统尚未建立需要了解未来系统的性能,只能通过对模型的试验来确认

  2. 在真实系统上进行试验可能会引起破坏或发生故障

  3. 系统无法恢复,如经济系统,若造成损失则无法挽回

  4. 试验条件无法保证,如进行多次试验难以保证试验条件相同,或试验时间太长,或费用昂贵

系统仿真的概念

仿真

系统仿真的必要性体现在以下方面:

  1. 优化设计:在复杂的系统建立以前,能够通过改变仿真模型结构和调整参数来优化系统设计,对系统和系统的某一部分进行性能评价
  2. 节省经费:仿真试验只需要在可重复使用的模型上进行,所花费的成本远比在实际产品上做经费低
  3. 故障诊断:系统发生故障后,设法使之重演,以便判断故障发生的原因
  4. 避免危险:某些试验有危险,不允许进行,而仿真试验可以避免危险性
  5. 假设预测:仿真可以预测系统的特性,也可以预测外部作用对系统的影响
  6. 可以训练系统操作人员
  7. 为管理决策和技术决策提供依据

系统仿真的分类

按仿真模型的种类分类
  1. 物理仿真:依据相似原理,把实际系统按比例缩小或放大制成物理模型,其状态变量和原系统完全相同
  2. 数学仿真:按照实际系统的数学关系构造系统的数学模型,并在计算机上进行实验研究
  3. 数学-物理仿真:将系统的物理模型和数学模型有机地组合在一起进行实验研究
按仿真模型与实际系统的时间关系分类
  1. 实时仿真:指仿真模型时钟τ\tau和实际系统时钟的比例关系为τt=1\displaystyle\frac{\tau}{t} = 1,是同步的,可以实时反映出实际系统的运行状态
  2. 超实时仿真:指仿真模型时钟τ\tau和实际系统时钟的比例关系为τt<1\displaystyle\frac{\tau}{t} < 1,即仿真模型的时钟超前于实际系统时钟
  3. 慢实时仿真:指仿真模型时钟τ\tau和实际系统时钟的比例关系为τt>1\displaystyle\frac{\tau}{t} > 1,即仿真模型的时钟滞后于实际系统时钟
按系统随时间变化的状态分类
  1. 连续系统仿真:系统的输入输出信号均为时间的连续函数,可用一组数学表达式来描述,此类系统可以用差分方程来描述
  2. 离散事件系统仿真:系统的状态变化只在离散时刻发生,且由某种随机事件驱动,多采用流程图或网络图表达

系统仿真的过程与特点

系统仿真的特征

系统仿真有三个基本要素:

  1. 实际系统
  2. 数学模型
  3. 计算机

三项基本活动:

  1. 模型建立
  2. 模型变换
  3. 仿真实验

系统仿真的过程

系统仿真的过程主要有以下几个方面:

  1. 系统定义:根据仿真目的,了解相关要求,确定所仿真系统的边界与约束条件等
  2. 数学建模:根据系统实验知识,仿真目的和实验数据来确定系统数学模型的框架,结构和参数,模型的繁简程度应与仿真目的相匹配,确保模型的有效性和仿真的经济性
  3. 仿真建模:根据数学模型的形式,计算机的类型以及仿真目的将数学模型变成适合于计算机处理的形式
  4. 模型输入:将仿真模型输入计算机,设定试验条件并进行记录
  5. 模型实验:根据仿真目的在模型上进行实验,即仿真
  6. 结果分析:根据试验要求对结果做分析,整理以及文档化,根据分析的结果修正数学模型,仿真模型,仿真程序,以进行新的仿真试验

系统仿真的特点

  1. 研究方法简单,方便,灵活,多样
  2. 实验成本低
  3. 实验结果充分

系统模型的建立与表示

模型建立

数学建模的概念

一个系统的数学模型可以定义为如下集合结构:

S=(T,X,Ω,Q,Y,δ,λ)S = (T,X,\Omega,Q,Y,\delta,\lambda)

其中:

TT表示时间基,描述变化的事件坐标,当TT为整数时称为离散事件系统,为实数时称为连续时间系统

XX表示输入集,代表外部环境对系统的作用

Ω\Omega表示输入段集,描述某个时间间隔内的输入模式,是(X,T)(X,T)的一个子集

QQ表示内部状态集,描述系统内部状态量,是系统内部结构建模的核心

YY表示输出集,系统通过它作用于环境

δ\delta表示状态转移函数,定义系统内部状态是如何变化的

λ\lambda表示输出函数,是一个映射

机理分析建模方法

一般是在若干简化假设条件下,以各学科专业知识为基础,对系统进行分析并确定模型结构,参数并列出相应的运动方程

使用机理建模法的前提是对系统的运行机理完全清楚,即白箱问题

机理分析建模法步骤如下:

  1. 分析系统功能,了解其结构和工作原理,对系统做出与建模目标相关的描述
  2. 找出系统的输入变量和输出变量
  3. 按照系统遵循的规律写出各部分的微分方程或传递函数等
  4. 消除中间变量,得到初步数学模型
  5. 进行模型标准化
  6. 进行模型验证

以上图的RLC电路为例,分析电路系统输入和输出之间的关系

输入为电源,输出为C上的电压

方法一:在时域中求解

根据KVL定律可得

e(t)=vR(t)+vL(t)+vo(t)e(t) = v_R(t) + v_L(t)+v_o(t)

再根据RLC元件的特点可得:

vR(t)=i(t)×Rv_R(t) = i(t)\times R

vL(t)=Ldi(t)dtv_L(t) = L\frac{di(t)}{dt}

i(t)=Cdvo(t)dti(t) = C\frac{dv_o(t)}{dt}

联立并整理上式可得:

LCd2vo(t)dt2+RCdvo(t)dt+vo(t)=e(t)LC\frac{d^2v_o(t)}{dt^2} + RC\frac{dv_o(t)}{dt} + v_o(t) = e(t)

输出变量在等号右侧,从左至右按阶次依次递减排列

方法二:利用拉氏变换进行求解

根据KVL定律可得

E(s)=VR(s)+VL(s)+Vo(s)E(s) = V_R(s) + V_L(s) + V_o(s)

再根据RLC元件的特点可得:

VR(s)=I(s)×RV_R(s) = I(s)\times R

VL(s)=SL×I(s)V_L(s) = SL\times I(s)

I(s)=sC×Vo(s)I(s) = sC\times V_o(s)

联立并整理可得:

s2LCVo(s)+sRCVo(s)+Vo(s)=E(s)s^2LCV_o(s) + sRCV_o(s) + V_o(s) = E(s)

H(s)=Vo(s)E(s)=1s2LC+sRC+1H(s) = \frac{V_o(s)}{E(s)} = \frac{1}{s^2LC + sRC + 1}

其中H(s)H(s)为系统的复频域模型

方法三:建立状态方程模型

在状态方程中,状态变量的选择不是特定的,可以任意选取(一般选择电感或电容上的电流或电压,作为状态变量)

本例中选择回路电流i(t)i(t)和输出电压vo(t)v_o(t)作为状态变量

根据KVL定律可得

e(t)=vR(t)+vL(t)+vo(t)e(t) = v_R(t) + v_L(t)+v_o(t)

再根据RLC元件的特点可得:

vR(t)=i(t)Rv_R(t) = i(t)R

vL(t)=Ldi(t)dt,di(t)dt=1LvL(t)v_L(t) = L\frac{di(t)}{dt},\frac{di(t)}{dt} = \frac{1}{L}v_L(t)

i(t)=Cdvo(t)dt,dvo(t)dt=1Ci(t)i(t) = C\frac{dv_o(t)}{dt},\frac{dv_o(t)}{dt} = \frac{1}{C}i(t)

整理可得:

i˙(t)=1Lvo(t)RLi(t)+1Le(t)\dot{i}(t) = -\frac{1}{L}v_o(t) - \frac{R}{L}i(t) + \frac{1}{L}e(t)

v˙o(t)=1Ci(t)\dot{v}_o(t) = \frac 1 Ci(t)

可以写成矩阵形式:

[v˙o(t)i˙(t)]=[01C1LRL][vo(t)i(t)]+[01L]e(t)\left [ \begin{array}{c} \dot v_o(t) \\ \dot i(t) \\ \end{array} \right ] = \left [ \begin{array}{cc} 0 & \frac 1 C \\ -\frac 1 L & -\frac R L \\ \end{array} \right ] \left [ \begin{array}{c} v_o(t) \\ i(t) \\ \end{array} \right ] + \left [ \begin{array}{c} 0 \\ \frac 1 L \\ \end{array} \right ] e(t)

上式即为该电路系统的状态方程模型

模型表示

系统的常微分方程模型

对于单输入单输出模型(SISO),其微分方程的一般形式为:

any(n)(t)+an1y(n1)(t)++a0y(t)=bmu(m)(t)+bm1u(m1)(t)++b0u(t)a_{n}y^{(n)}(t) + a_{n - 1}y^{(n - 1)}(t) + \cdots + a_0y(t) = b_{m}u^{(m)}(t) + b_{m - 1}u^{(m - 1)}(t) + \cdots + b_0u(t)

其中y,uy,u为系统的输出和输入,ai,bja_i,b_j表示各项输出和输入的系数,均为常系数

num=[bm,bm1,,b1,b0]num = [b_m,b_{m - 1},\cdots,b_1,b_0]表示方程等号右侧的各项系数

den=[an,an1,,a1,a0]den = [a_n,a_{n - 1},\cdots,a_1,a_0]表示方程等号左侧的各项系数

缺项均用0补齐

若系统的输入和输出不是一个,而是多个,则称为多输入多输出系统(MIMO)

系统的传递函数模型

传递函数只与系统本身的结构,特性和参数有关,而与输入输出量的变化无关

对于SISO连续系统,系统的传递函数如下所示:

G(s)=B(s)A(s)=bmsm+bm1sm1++b1s+b0ansn+an1sn1++a1s+a0(nm)G(s) = \frac{B(s)}{A(s)} = \frac{b_ms^m+b_{m - 1}s^{m - 1} + \cdots + b_1s+b_0}{a_ns^n + a_{n - 1}s^{n - 1} + \cdots + a_1s + a_0}\quad(n\ge m)

在MATLAB当中,用函数tf可以建立一个连续系统传递函数模型,其调用格式为:

1
sys = tf(num,den)

其中,num为传递函数分子系数向量,den为传递函数分母系数向量

对于SISO离散时间系统,对传递函数进行z变换,则可得到该离散系统的脉冲传递函数:

G(z)=Y(z)U(z)=fmzm+fm1zm1++f0gnzn+gn1zn1++g0G(z) = \frac{Y(z)}{U(z)} = \frac{f_mz^m + f_{m - 1}z^{m - 1} + \cdots + f_0}{g_nz^n + g_{n - 1}z^{n - 1} + \cdots + g_0}

其中,对于LTI离散系统而言,fif_igig_i均为常数

在MATLAB中可以用tf函数来建立系统的函数模型,调用格式为:

1
sys = tf(num,den,Ts)

其中num为z传递函数分子系数向量,den为z传递函数分母系数向量,Ts为采样周期

在MATLAB中也可以用多项式乘法函数conv来描述传递函数,例如

G(s)=4(s+2)(s2+6s+6)2s(s+1)3(s3+3s2+2s+5)G(s) = \frac{4(s + 2)(s^2 + 6s + 6)^2}{s(s + 1)^3(s^3 + 3s^2 + 2s + 5)}

可以用如下语句进行描述:

1
2
3
num = 4 * conv([1,2],conv([1,6,6],[1,6,6]));
den = conv([1,0],conv([1,1],conv([1,1],conv([1,1],[1,3,2,5]))));
sys = tf(num,den)

系统的状态空间模型

微分方程和传递函数均是描述系统性能的数学模型,它只能反映出系统输入和输出之间的对应关系,通常称之为外部模型,而在系统仿真时,常常要了解系统中各内部变量的状态,这就需要用到状态空间模型

在状态空间模型中,状态变量既可以是系统的输入和输出变量,也可以是系统的内部变量

给定一个系统,当引入nn个状态变量时,可得到由nn个一阶微分方程组成的方程组,采用矩阵描述,可以得到:

x˙=Ax+By\pmb{\dot x} = \pmb{Ax} + \pmb{By}

y=Cx+Du\pmb y = \pmb{Cx} + \pmb{Du}

其中,x\pmb x为状态向量,u\pmb u为输入向量,y\pmb y为输出向量

A\pmb A为状态变量系数矩阵,称为系统矩阵

B\pmb B为输入变量系数矩阵,称为输入矩阵

C\pmb C为输出变量系数矩阵,称为输出矩阵

D\pmb D称为直接传递矩阵

x˙=Ax+By\pmb{\dot{x}} = \pmb {Ax} + \pmb {By}为系统的状态方程,y=Cx+Du\pmb y = \pmb{Cx} + \pmb{Du}为系统的输出方程,两者组合后为系统的状态方程

在MATLAB中,用ss函数可以建立一个连续系统状态空间模型,调用格式为:

1
sys = ss(A,B,C,D)

其中,ABCD,为系统状态方程系数矩阵

对于离散时间系统而言,状态空间模型可以写成:

X(k+1)=FX(k)+GU(k)\pmb X(k + 1) = \pmb{FX}(k) + \pmb{GU}(k)

Y(k+1)=CX(k+1)+DU(k+1)\pmb Y(k + 1) = \pmb{CX}(k + 1) + \pmb{DU}(k + 1)

在MATLAB中,用ss函数也可以建立一个离散时间系统的传递函数模型,其调用格式为:

1
sys = ss(A,B,C,D,Ts)

其中,ABCD,为系统状态方程系数矩阵,Ts为采样周期

系统模型的其他形式

零、极点增益模型(Zero-Pole)

零、极点增益模型是通过对系统传递函数分子和分母多项式进行分解,求得系统的零、极点而得到,即将传递函数用常数项与系统的零、极点来表示,这个常数项通常记作kk,称作系统的增益

对于SISO连续系统而言,其零、极点模型为:

G(s)=k(sz1)(sz2)(szm)(sp1)(sp2)(spn)G(s) = k\frac{(s - z_1)(s - z_2)\cdots(s - z_m)}{(s - p_1)(s - p_2)\cdots(s - p_n)}

其中,ziz_ipjp_j为系统的零点和极点,kk为系统增益

在MATLAB中,可以用函数zpk来直接建立连续系统的零、极点增益模型,其调用格式为:

1
sys = zpk(z,p,k)

其中,z,p,k分别表示系统的零点向量,极点向量和增益

对于离散时间系统,也可以用函数zpk建立零、极点模型,其调用格式为:

1
sys = zpk(z,p,k,Ts)

其中Ts为采样周期

对于复杂的,不便于求解零、极点的系统,可以用MATLAB的求根函数root,调用格式为:

1
2
z = root(num)
p = root(den)

其中numden分别为传递函数模型的分子和分母多项式系数向量

对于MIMO系统,函数zpk也可以建立其零、极点增益模型,调用格式和SISO相同,但是zpk不再是一维向量,而是矩阵

频率响应数据模型(Frequency-Response-Data)

当我们无法直接建立系统的传递函数或状态空间模型,而只知道该系统在某些频率点上的响应时,可以利用函数frd建立该系统的频率响应模型,其调用格式为:

1
G = frd(resp,freq)

其中,系统的频率响应是复数,可用response=[g1,g2,,gk]response = [g_1,g_2,\cdots,g_k]输入,对应的频率ω\omegafreq=[ω1,ω2,,ωk]freq = [\omega_1,\omega_2,\cdots,\omega_k]输入,两者应有相同的列数

模型转换

系统的模型转换

各种数学模型都有其特点和适用的场合和领域,所以有时候需要进行模型的转换

LTI的模型有TF,ZPK,SS和FRD4类,每类模型又有连续和离散模型之分

函数 说明
residue 由传递函数形式转化为部分分式形式
ss2tf 由状态空间形式转化为传递函数形式
ss2zp 由状态空间形式转化为零、极点形式
tf2ss 由传递函数形式转化为状态空间形式
tf2zp 由传递函数形式转化为零、极点形式
zp2ss 由零、极点形式转化为状态空间形式
zp2tf 由零、极点形式转化为传递函数形式
c2d 将状态空间模型由连续形式转化为离散形式
c2dm 连续形式到离散形式的转换
c2dt 连续形式到离散形式的对输入纯时间延迟转换
d2c 将状态空间模型由离散到连续的转换
d2cm 离散形式到连续形式的方法转换

TF,ZPK和SS模型之间可以相互转换

FRD模型不可以直接转换成TF,ZPK和SS模型,但是可以由其他3种模型转换得到

系统模型的连接

模型串联

函数series用于两个线性模型串联,其调用格式为:

1
sys = series(sys1,sys2)

其中,sys1sys2sys如图所示

对于SISO模型,series函数会将第一个子系统的输出连接到第二个子系统的输入,其用法为:

1
[A,B,C,D] = series(A1,B1,C1,D1,A2,B2,C2,D2)

对于MIMO模型,函数series的调用格式为:

1
sys = series(sys1,sys2,outputs1,inputs2)

该函数在执行连接sys1sys2时会将第一个子系统的第一个输出接到第二个子系统的输入上

系统端口的名称可以用函数set设置

暂略

系统的仿真分析

系统的稳定性分析

稳定是系统能够正常运行的首要条件

稳定性的一般定义是:设一线性定常系统原处于某一平衡状态,若它瞬间受到某一扰动作用而偏离了原来的平衡状态,当此扰动撤销后,系统仍能回到原有的平衡状态,则称该系统是稳定的,反之则为不稳定

线性系统的稳定性取决于系统的固有特征,而与系统的输入信号无关

间接判别法

线性系统稳定的充要条件是闭环特征方程式的根必须都位于ss的左半平面

令系统的闭环特征方程为:

α0sn+α1sn1+α2sn2++αn1s+αn=0α0>0\alpha_0s^n + \alpha_1s^{n - 1} + \alpha_2s^{n - 2} + \cdots + \alpha_{n - 1}s + \alpha_n = 0\quad\alpha_0 > 0

将各项系数排成劳斯表:

sna0a2a4a6sn1a1a3a5a5sn2b1b2b3b4sn3c1c2c3s2d1d2d3s1e1e2s0f1\begin{align} s^n& \qquad a_0 \quad a_2 \quad a_4 \quad a_6 \quad \cdots\\ s^{n - 1}& \qquad a_1 \quad a_3 \quad a_5 \quad a_5 \quad \cdots\\ s^{n - 2}& \qquad b_1 \quad b_2 \quad b_3 \quad b_4 \quad \cdots\\ s^{n - 3}& \qquad c_1 \quad c_2 \quad c_3 \quad \cdots\\ \vdots\\ s^2& \qquad d_1 \quad d_2 \quad d_3\\ s^1& \qquad e_1 \quad e_2 \\ s^0& \qquad f_1 \end{align}

其中

b1=a1a2a0a3a1b2=a1a4a0a5a1b3=a1a6a0a7a1c1=b1a3a1b2b1c2=b1a5a1b3b1c3=b1a7a1b4b1f1=e1d2d1e2e1\begin{align} b_1& = \frac{a_1a_2-a_0a_3}{a_1}\quad b_2 = \frac{a_1a_4-a_0a_5}{a_1} \quad b_3 = \frac{a_1a_6-a_0a_7}{a_1}\quad\cdots\\ c_1& = \frac{b_1a_3-a_1b_2}{b_1} \quad c_2 = \frac{b_1a_5-a_1b_3}{b_1} \quad c_3 = \frac{b_1a_7-a_1b_4}{b_1}\quad \cdots\\ \vdots&\\ f_1& = \frac{e_1d_2-d_1e_2}{e_1} \end{align}

凡在运算中出现的空位用0填充,一直计算下去,就可以得到n+1n+1行系数

判决规则如下:

  1. 如果所列劳斯表第一列系数均为正值,则其特征方程式的根都在SS域左平面,相应系统是稳定的
  2. 如果所列劳斯表第一列系数的符号有变化,其变化的次数等于该特征方程式的根在SS域右半平面的个数,相应的系统是不稳定的

也可以通过直接判别法求稳定性,即设系统的传递函数的分母为0,求解所有零极点

若所有根在左半平面/单位圆内(拉氏变换/zz变换),则系统稳定

自动控制系统概述

自动控制的概念

自动控制是指在没有人直接参与的情况下,利用外加的设备或装置(控制装置),使机器、设备或生产过程(控制对象)的某个工作状态或参数(被控量) 自动地按照预定的规律运行

自动控制的方式

  1. 按给定值操纵的开环控制

开环控制:系统的输入端和输出端之间不存在反馈回路,输出量对系统的控制作用没有影响

  1. 按干扰补偿的开环控制

利用干扰信号产生控制作用,以及时补偿干扰对被控量的直接影响

  1. 按偏差调节的闭环控制

通过对系统的输出量进行测量,利用实际值和误差值的差去对系统进行调整

  1. 复合控制

PID控制原理

PID控制概述

输入:控制偏差e(t)=r(t)y(t)e(t) = r(t) - y(t)

输出:偏差的比例§、积分(I)和微分(D)的线性组合

u(t)=Kc(e(t)+1TI0te(t)dt+TDde(t)dt)u(t) = K_c(e(t) + \frac{1}{T_I}\int^t_0e(t)dt + T_D\frac{de(t)}{dt})

其中KcK_c为比例系数,TIT_I为积分时间常数,TDT_D为微分时间常数

基本的PID控制规律可以描述为:

Gc(s)=Kp+KIs+KDsG_c(s) = K_p + \frac{K_I}{s} + K_Ds

控制方式介绍

比例§控制

输入输出成比例关系

Gc(s)=KpG_c(s) = K_p

其中KpK_p为比例系数或增益

在过程控制中,通常用比例度δ\delta表示控制输出与偏差成线性关系的比例控制器输入(偏差)的范围

δ=1Kc×100%\delta = \frac{1}{K_c}\times 100\%

比例微分(PD)控制

传递函数为:

Gc(s)=Kp+KpτsG_c(s) = K_p+K_p\tau s

其中KpK_p为比例系数,τ\tau为微分时间常数

输出信号为:

u(t)=Kpe(t)+Kpτde(t)dtu(t) = K_pe(t) + K_p\tau\frac{de(t)}{dt}

比例积分(PI)控制

传递函数为:

Gc(s)=Kp+KpTi1s=Kp(s+1Ti)sG_c(s) = K_p + \frac{K_p}{T_i}\cdot\frac{1}{s} = \frac{K_p(s + \frac{1}{T_i})}{s}

其中KpK_p为比例系数,TiT_i为积分时间常数

输出信号为:

u(t)=Kpe(t)+KpTi0te(t)dtu(t) = K_pe(t) + \frac{K_p}{T_i}\int^t_0e(t)dt

比例积分微分(PID)控制

输出信号为:

u(t)=Kce+S00tedt+S2dedtu(t) = K_ce + S_0\int^t_0edt + S_2\frac{de}{dt}

其中的三个特征参数为:比例带δ\delta、积分时间TiT_i、微分时间TdT_d

PID参数整定方法
  1. Ziegler-Nichols整定方法

假设对象的近似传递函数为:G(s)=KeLsTs+1\displaystyle G(s) = \frac{Ke^{-Ls}}{Ts + 1}

其中LL为延迟时间,KK为放大系数,TT为时间常数

可得参数如下表

控制器类型 比例度δ\delta 积分时间TiT_i 微分时间τ\tau
P T/(KL)=0.5KpsT / (KL) = 0.5K_{ps} \infin 00
PI 0.9T/(KL)=0.45Kps0.9T/(KL) = 0.45K_{ps} L/0.3L/0.3 00
PID 1.2T/(KL)=0.6Kps1.2T/(KL) = 0.6K_{ps} 2.2L2.2L 0.5L0.5L
  1. 临界比例度法

设临界比例度为δk\delta_k,临界振荡周期为TkT_k

采用临界比例度法的前提是系统的阶数在3阶或3阶以上

控制器类型 比例度δ\delta 积分时间TiT_i 微分时间τ\tau
P 2δk2\delta_k \infin 00
PI 2.2δk2.2\delta_k 0.833Tk0.833T_k 00
PID 1.7δk1.7\delta_k 0.50Tk0.50T_k 0.125Tk0.125T_k
  1. 衰减曲线法