混合线性模型（精选9篇）
混合线性模型 篇1
水稻在我国有着广泛的种植面积, 稻田持续健康增产对我国的经济发展起重要作用。水稻产量受多方面因素影响, 如气候, 地形, 品种和土壤肥力等[1], 其中土壤肥力, 表征着土壤的水、肥、气、热状况, 直接影响作物的生长。研究表明, 合理施肥和管理不仅能够提高作物产量, 而且可以减少环境损害[2], 而土壤肥力是土壤地力的基础, 高肥力土壤表明高的土地生产力, 即高地力。
土壤地力对产量的贡献不容忽视, 唐旭等[3]得到水稻大麦轮作体系下土壤地力对大麦和水稻产量的贡献率分别为69%和75%-81%, 马常宝等[4]认为潮土区土壤地力对小麦和玉米产量的贡献分别达51.4%和54%, 国外研究也表明高地力土壤能在一定时期内维持产量的持续稳定增长[5]。因此, 对土壤地力合理评价有助于了解土地的生产力, 从而采取相应措施提升土壤地力以尽可能发挥土壤的增产潜能。土壤地力评价比较容易的是通过施肥和不施肥产量对比或依据地力调查来评估[6,7], 目前较多的是采用更为复杂的模型方法进行评价[8,9,10], 骆东奇等[11]对土壤肥力评价指标和方法进行了总结, 但基于剖面上土壤理化性质的多元混合线性回归法则不多见, 该法考虑到不同层次上理化性质间的差异性和相互作用, 引入随机效应使得预测将更加精准[12]。
世界上最早的水稻种植区发现于湖南澧县[13], 而在道县则发现了迄今为止最早的古栽培稻谷[14], 因此, 本研究在两地分别采集了高中低三种产量的稻田土壤剖面样品, 分层对其理化性质进行测定, 然后建立稻田地力与理化性质间的多元线性回归模型和多元混合线性模型并对其进行比较分析, 同时比较不同土壤深度层次划分采样时所建立的多元混合线性模型优劣, 用于更好指导取样分析, 从而期望为稻田地力评估和指导农业生产实践提供理论参考。
1 材料与方法
1.1 研究区概况
本文选取的稻田样点分别位于湖南澧县和道县。澧县属亚热带湿润季风气候, 平均海拔59 m。该县春温秋和, 夏热冬寒, 雨量充足, 年均气温16.5℃, 年降水量1213 mm, 日照年平均1770 h, 无霜期265 d。全境以水稻土, 红壤和潮土为主, 地貌特征主要为平原和岗地, 农业种植以水稻、棉花、油料作物为主。道县也属亚热带湿润季风气候区, 平均海拔181 m。境内热量丰富, 无霜期长, 雨量充沛, 光照充足, 光能潜力大。境内年均气温18.6℃, 无霜期302 d, 年降雨量1523 mm, 年日照时数1514 h, 年蒸发量1586 mm。境内季风气候明显, 冬季盛吹北风, 气温较低;夏季盛吹偏南风, 气温较高;春秋两季温度适中, 因此四季分明。
1.2 样品采集与测定
利用土壤柱状圆筒采样器 (Eijkelkamp, 荷兰, 内径9 cm) 和剖面挖掘方法于2011年11月分别采集了两县稻田的高中低产50 cm深土壤剖面样品各一个。澧县高产田为稻-油轮作, 水稻品种为徽两优3号, 油菜品种为沛油520, 一季稻高中低产田稻谷产量分别为10.65、7.5、5.25 t/hm2。道县稻田种植两季, 高中低产田单季产量分别为7.5、6.0、5.25 t/hm2。
对于每一个土壤剖面, 按从上到下每隔5 cm取原状土环刀 (100 cm3) 并相应取该层次土壤剩余样品, 环刀样用于测定土壤两个物理性质, 即土壤饱和导水率 (Ks) 和容重 (BD) , 对应层次的剩余土样风干后, 分别过2 mm和0.25 mm筛用于测定土壤理化性质, 包括物理性质:机械组成;化学性质:阳离子交换量 (CEC) 及单个的钾 (K+) 、钙 (Ca2+) 、钠 (Na+) 和镁 (Mg2+) 离子浓度, p H, 电导率 (EC) ;养分特征:总氮 (TSN) , 总磷 (TSP) , 总钾 (TSK) , 有机碳 (SOC) , 速效磷 (AVP) , 速效钾 (AVK) 和碱解氮 (AVN) 。
CEC采用氯化钡提取-快速滴定法[15], AVN采用碱解扩散法[16], SOC利用碳氮分析仪测定, Ks采用双环刀恒定水头法, 土壤理化性质测定主要参照刘光崧[17]主编的《土壤理化分析与剖面描述》, 机械组成利用比重计法测定, 采用美国制, 即砂粒 (SAND) :0.25-2.0 mm;粉粒 (SILT) :0.25-0.002 mm;粘粒 (CLAY) <0.002 mm。
1.3 数据处理与统计
1.3.1 数据处理
本文依据所选稻田水稻产量情况, 将稻田地力定量为不同数值, 即高产=3;中产=2;低产=1。若认为土壤理化性质在同一土壤剖面中不同深度层次上相互独立, 将所有理化性质与稻田地力作相关性分析, 选取显著的影响因子, 利用逐步回归法构建多元线性回归模型。
若考虑在同一剖面中, 不同深度层次同一土壤理化性质间具有关联性, 则对土壤深度分成四组进行考虑, 分别为:10层, 即0-50 cm按5 cm进行分层;5层:即0-50 cm按10 cm进行分层;3层:即为0-20 cm, 20-30 cm, 30-50 cm;分2层, 即按犁底层划分为0-20 cm和20-50 cm。由于分层情形不同, 每层中所含的实测样品数则不同, 若分10层, 则每层含有6个样品数据;分5层, 则每层含有12个数据;分3层, 则对应层样品数分别为24, 12, 24;分2层, 则对应层样品数分别为24和36。利用每层中样本数据进行相关性分析, 选出在每层上都具有显著性的土壤理化性质, 然后分别建立四个多元混合线性回归模型并比较, 从而得到预测效果最好的分层采样方法用于指导取样分析。
1.3.2 统计方法
依据不同假设, 本研究分别建立多元线性回归模型和Laird和Ware[12]提出的多元混合线性模型来预测稻田地力并比较。
多元线性回归模型一般表达式:
式中:Y为因变量向量, X为自变量向量, β为回归系数向量, ε为残差向量。
多元混合线性模型一般表达式:
式中:i为分组, Yi表示第i组因变量向量, Xi为第i组固定效应向量, βi为第i组固定效应系数向量, Zi为第i组随机效应向量, λi为第i组随机效应系数向量, εi为第i组残差向量。
多元混合线性模型, 其实质是对多元线性回归模型的残差项作进一步的拆解, 进而引出随机效应, 它能够获得更准确的预测效果。多元混合线性模型常用于对聚类或非独立的分层数据以及同一个体不同时间段的重复测量数据进行拟合[18], 该模型同时纳入了对不同聚类或分层数据产生的固定效应和随机效应影响, 固定效应都是相同的, 不同只是随机效应。随机效应主要源于从总体中取样时带来的波动性变化, 对样品数据的影响具有非系统性和难以预测性, 而固定效应则相反, 对样品数据具有系统性和可预见性特征[19]。
利用均方根误差 (RMSE) 和调整的确定系数 (R2adj) 来评估多元线性回归模型和多元混合线性模型预测优劣, RMSE值越小, 表明模型预测越准确。R2adj位于0到1之间, 其值越接近1, 模型拟合越准。而针对不同分组构建的不同多元混合线性回归模型比较则采用方差分析, 其评判依据是赤池信息准则 (Akaike information criterion, AIC) , 贝叶斯信息准则 (Bayesian information criterion, BIC) 和残差 (ε) , 这些值越小, 表明模型拟合越好。评判指标的计算公式分别表示为:
式中:ym和yp分别为样品观测值和预测值, n为观测样品数。R2为模型确定系数, M为变量个数。k为模型中未知参数个数, L为模型最大似然函数值,
利用R统计软件 (http://www.r-project.org/) 进行相关性分析、作图和模型构建。多元混合线性模型中选用的软件包为lmer 4.0, 建模前进行相关性分析并分析不同分组情形下相关系数分布, 然后提取显著的固定效应 (P<0.05) 和随机效应进行建模, 不同分组情况下查表得到的对应相关系数临界值分别为:10层:0.755 (n=6) ;5层:0.553 (n=12) ;3层:0.553 (n=12) ;2层:0.396 (n=24) , 即认为相关系数大于该临界值时表明相关性显著, 多元混合线性模型构建过程中会自动计算出相应的AIC、BIC和ε。
考虑到Ks具有较强空间变异性, 通常符合对数正态分布[20], 因此本文对Ks进行了常用对数转换, 即K10=log10 (Ks) 。
2 结果与分析
2.1 土壤性质
高中低产稻田土壤剖面样品所有理化性质特征见表1, 土样质地以粉壤土和壤土为主, 分别占45%和32%。在三种不同地力稻田土壤中, K10都表现出较强的空间变异性, 且在高产田中变异系数最大, 但K10的平均值在高产田中最小, AVP在高产田中也表现为强变异 (CV>100) , 而其他性质在不同地力土壤中都属于中等变异 (10<CV≤100) 。高产田中的SILT, CLAY, SOC, AVK, p H, CEC, Ca2+和K+较中、低产田都相对高, 而BD, TSP, TSK, AVN, AVP, EC, Mg和BCSR在中产田中平均值较高, 只有SAND在低产田中最大, TSN含量在高中低产田中差异不明显。
注:K10:饱和导水率;BD:容重;CLAY/SILT/SAND:粘粒/粉粒/砂粒含量;SOC:有机碳;TSN:总氮;TSP/TSK:总磷/总钾;AVN/AVP/AVK:碱解氮/速效磷/速效钾;p H:p H值;Ca2+/K+/Mg2+/Na+/CEC:钙/钾/镁/钠/阳离子交换量;EC:电导率;BCSR:盐基饱和度;CV:变异系数 (%) 。
2.2 因子筛选
从多元线性回归模型中的土壤理化性质与稻田地力间相关系数可知 (表2) , 影响地力的显著土壤性质中, Ca2+和CEC影响最大, 相关系数都为0.59 (P<0.01) , 其次为K+和AVK, 相关系数分别为0.44 (P<0.01) 和0.39 (P<0.01) , 最后为p H和CLAY, 其他性质都不显著。Na+, EC, BCSR, K10和SAND表现出与地力的负相关性, 其他性质都为正相关性。
注:YIELD_N:稻田地力;K10:饱和导水率 (cm/d) ;BD:容重 (g/cm3) ;CLAY/SILT/SAND:粘粒/粉粒/砂粒含量 (%) ;SOC:有机碳 (%) ;TSN:总氮 (%) ;TSP/TSK:总磷/总钾 (g/kg) ;AVN/AVP/AVK:碱解氮/速效磷/速效钾 (mg/kg) ;p H:p H值;Ca2+/K+/Mg2+/Na+/CEC:钙/钾/镁/钠/阳离子交换量 (cmol/kg) ;EC:电导率 (u S/cm) ;BCSR:盐基饱和度 (-) 。*:P<0.05;**:P<0.01。
从多元混合线性模型中的土壤理化性质与稻田地力间相关系数随土壤深度变化可知 (图1) , 化学性质方面, Ca2+, CEC和p H表现为从上层至下层的同步增长趋势且都显著, 并且Ca2+和CEC在20-30cm处呈现出于地力最强的相关性, 而K+则呈现同步下降趋势且都显著, Mg2+, EC和BCSR在不同分层上对稻田地力影响都不显著, Na+在土壤深度10-30 cm处表现出与稻田地力较强的负相关性。养分特征中只有AVP, AVK和TSP在不同层次上与地力关系都紧密 (P<0.05) , 从土壤剖面上层到下层, AVK都呈现正相关, 分布比较均匀, 对于TSP则在上层相关性比在下层强, 而AVP在20-30 cm处与地力相关显著性不明显。AVN、TSN和TSK与地力的相关系数在各深度层次上都较低且在不同层次上的相关性呈正负交替变化, SOC在20-30 cm处呈现较强的正相关性, 再往下则变为负相关。物理性质中, K10在上层与稻田地力为显著正相关, 20 cm以下则呈现显著负相关。BD与K10有着较为相似的相关系数分布趋势。CLAY在下层对稻田地力影响较大且为正相关, 而SAND则在下层表现出与地力较强的负相关。通过对图1的分析, 相关性从上层至下层转变剧烈且对稻田地力影响较大的性质有K10, SOC, Na+和BD, 这些性质在多元线性回归建模中都没有表现出较好的相关性, 但是经过分层处理后则在不同层次上会有显著性变化出现, 表明这些性质受分层影响, 建立混合效应模型时将作为随机效应, 而其他显著变量如Ca2+, CEC, p H, K+, AVP, AVK, TSP, CLAY和SAND则将作为固定效应进行模型构建。
2.3 模型分析
建立的多元线性回归模型如下:
该模型只选用了四个变量, 分别为p H、K+、Ca2+和TSP, 能解释稻田地力总变异的73%, 这说明还有一些未知的因素未考虑到而影响预测精度, 方程中K+和TSP回归系数较高, 表明模型对其变化比较敏感。模型中未选用显著变量CEC和AVK, 而选用p H和TSP, 主要是由于Ca2+和CEC, 以及K+和AVK之间的显著相关性, 在多元逐步回归计算中进行了筛选, 同时考虑到其他较为重要的性质。
不同分组情形下建立的多元混合线性模型对稻田地力的预测比较得出, 土壤深度分为3层, 即将土壤剖面按0-20、20-30、30-50 cm分层采样所构建的多元混合线性模型预测效果最好, 此时AIC和BIC值都最小 (表1) , R2adj达到0.85, 表示该模型能解释稻田地力的85%。与多元线性回归比较时, 多元混合线性模型的RMSE较小且R2adj较高, 表明多元混合线性模型预测比多元线性回归预测更准确。
土壤分3层构建的稻田地力多元混合线性模型包括固定效应 (表4) 和随机效应 (表5) , 选取的固定效应包括AVP、AVK、TSP、Ca2+、K+、CEC、p H、CLAY和SAND;随机效应包括K10、SOC、Na+和BD。总变异中由Na+和BD所带来的变异比值较大, 分别占模型总变异值的64.5%和9.7%, 表明这两个变量对稻田地力的多元混合线性模型会产生更大的随机影响。
注:Ca2+/K+/Mg2+/Na+/CEC:钙/钾/镁/钠/阳离子交换量 (cmol/kg) ;p H:p H值 (-) ;EC:电导率 (u S/cm) ;BCSR:盐基饱和度 (-) ;AVN/AVP/AVK:碱解氮/速效磷/速效钾 (mg/kg) ;TSN:总氮 (%) ;TSP/TSK:总磷/总钾 (g/kg) ;SOC:有机碳 (%) ;K10:饱和导水率 (cm/d) ;BD:容重 (g/cm3) ;CLAY/SILT/SAND:粘粒/粉粒/砂粒含量 (%) ;虚线为对应相关性显著临界值:即10层:0.755 (n=6) ;5层:0.553 (n=12) ;3层:0.553 (n=12) ;2层:0.396 (n=24) 。
利用所建立的多元混合线性模型和多元线性回归模型观测值与预测值作箱线图比较分析 (图2) , 多元混合线性模型对稻田地力预测值相对多元线性回归模型都较为集中, 尤其是对中、低地力的预测更好, 而在高地力预测方面, 简单模型预测值虽然表现较为集中, 但存在较多异常值, 进而影响了模型整体估测精度。
3 讨论
3.1 土壤理化性质对稻田地力的影响
土壤地力表征着土壤的基础潜能, 反映了土地的出产力。土壤肥力是土壤地力的重要属性, 高肥力土壤意味着高地力。在高地力土壤上往往在相对较少的人为管理或施肥下都能获得较高作物产量[3]。本研究将高中低产稻田的地力情况定量化, 即高产田代表高地力土壤, 中产田代表中等地力, 低产田代表低地力土壤。许多研究者认为土壤自身理化性质及肥力特性对作物产量起重要作用[21,22], 本研究发现高产田中土壤平均SILT, CLAY, SOC, AVK, p H, CEC和Ca2+含量较高, 而K10最低。SILT和CLAY代表土壤颗粒大小和含量, 细质地土壤有助于保持土壤养分进而提升土壤肥力。从图1中看出, K10在表层 (0-20 cm) 表现为与地力较强的正相关性, 底层为负显著相关性。表层K10越大, 土壤导水性越强, 说明土壤空隙较充分且土壤较为疏松, 这有助于养分运输和分配, 而底层若K10越大, 则土壤将难以保持充足养分和水分, 表明高地力土壤应在表层有较高的导水性而在底层则导水率应小。
注:Int:截距;AVP/AVK:速效磷/速效钾 (mg/kg) ;TSP:总磷 (g/kg) ;Ca2+/K+/CEC:钙/钾/阳离子交换量 (cmol/kg) ;p H:p H值 (-) ;CLAY/SAND:粘粒/砂粒含量 (%) 。
注:Int:截距;K10:饱和导水率 (cm/d) ;SOC:有机碳 (%) ;Na+:钠离子 (cmol/kg) ;BD:容重 (g/cm3) 。
图2多元线性回归模型和多元混合线性回归模型对稻田地力的预测值与观测值的比较Fig.2 Comparison between the observed and predicted soil productivity by the multiple linear regression model (left panel) and the mixed-effects linear regression model (right panel)
养分性质方面, SOC表征着土壤有机质, 吴槐泓等[23]通过定位试验认为有机质能指示土壤肥力水平并与水稻产量成正相关。本文从多元线性回归模型中的相关性分析发现SOC与地力的相关性很弱 (表2) , 但经过分层处理后, 发现在土壤20-30 cm深度处SOC与地力表现较强的正相关性。SOC在土壤中的含量一般随深度而逐渐降低, 但本研究认为犁底层处较高的SOC含量有助于提升稻田地力。AVK在高产田中含量高说明提高土壤中有效钾含量有助于提升地力 (图1) , 蒋德安等[24]认为适当供钾可以显著提高水稻谷粒产量, 这表明土壤中钾素的不足将限制着土壤地力, 但从图1发现, 土壤TSK并不与地力呈现显著相关, 只有土壤中AVK和K+增加, 才能显著提升稻田地力。
化学性质方面, 土壤深度10 cm以下, p H与地力都呈现显著正相关 (图1) , 表明底层高p H有助于增强稻田地力, 杨曾平[25]通过室内盆栽试验认为以石灰岩风化物发育的水稻土肥力比河流沉积物和第四纪红色粘土发育的水稻土高。CEC和Ca2+含量在高产田中表现比较突出, 与稻田地力呈显著正相关, 表明高地力稻田对阳离子的需求量比较大, 尤其是钙离子。本研究发现低产田中砂粒含量较大, 这跟砂粒不利于养分和水分的保持有关。通过对比高中低产稻田土壤理化性质分析发现, 稻田地力不受土壤总氮的影响, 这可能与目前的施肥制度有关, 通过施肥提高产量已成为目前农田中最常用的手段, 因此不同水稻产量稻田中总氮量已趋于饱和, 可能实际上对于高地力土壤本不需要施入当前用量的氮肥, 多余的氮则将导致氮肥浪费和环境污染。研究还发现某些土壤性质在考虑分层影响后, 与地力的相关性在某些层次上呈现显著相关性, 如Na+在20-30cm处表现较强的负相关性, K10, BD和SOC表现为表层较强的正相关而底层较强的负相关性 (图1) , 表明这些土壤性质对地力的影响受分层影响。
3.2 稻田地力模型
多元线性回归得到的稻田地力预测模型只包含了p H, K+和Ca2+以及TSP四个变量, 变量相对较少, 但模型预测精度只有73%, 并且多元线性回归模型要求样本数据需具有独立性, 然而考虑到土壤剖面样品在不同层次间可能存在相互影响, 如SOC随深度减小变化, 上层SOC与下层SOC可能具有较好的相关性, 因此建立了预测稻田地力的多元混合线性模型, 即同时选取固定效应和随机效应, 固定效应选取与地力具有显著相关性但基本不随土壤分层而发生正负相关性变化的土壤性质, 而随机效应则选那些因分层而导致的显著性及相关性正负变化的性质, 随机效应在不同层次上对地力总体预测表现并不一致, 由于混合效应考虑到了更多变异来源, 对重复数据或聚类或分层数据具有较好的应用前景[26,27], 理论上多元混合线性模型将提高预测精度。本文得到的多元混合线性模型将预测精度从多元线性回归模型的73%提高到了85%, 但多元混合线性模型选用了比较多的固定效应, 这将增加实验分析难度。因此, 应对于某些变量进一步筛选, 如TSP, K+, AVP, SAND和CLAY对地力的影响在不同层次影响其实并不都完全显著, 同时, 可能还有一些更好的变量未考虑, 如王远敏[28]通过田间控制试验认为稻田施硅肥有利于产量增加。
4 结论
本研究发现不同地力稻田土壤总氮差异不明显, 低产田中砂粒平均含量偏高而其他土壤理化性质都偏低。通过混合效应模型中相关性分析筛选出13个预测稻田地力的因素, 即固定效应9个, 包括速效钾、速效磷、总磷、钙离子、钾离子、阳离子交换量、p H值、粘粒和砂粒含量, 其中阳离子交换量与钙离子浓度对稻田地力影响最显著且正相关。随机效应4个, 即饱和导水率, 有机碳, 钠离子和容重, 其中钠离子浓度产生的随机效应最强。建立的多元线性回归模型对钾离子浓度和总磷含量变化比较敏感, 本文建立的稻田地力混合效应模型较多元线性回归模型预测性能从73%提高到了85%。
由于研究无法充分考虑到其他可能对稻田地力有更大影响的因素, 预测依旧存在一定的偏差, 后续可通过增加样品量以及选取更好的自变量以提高多元混合线性模型的预测精度和实用性能。
混合线性模型 篇2
基于Lyapunov第二方法,讨论标称中立型时滞系统的H∞性能问题.在此基础上,研究了具混合变时滞以及不确定性具时变结构的.中立型系统的鲁棒H∞控制器的设计方法, 推广和改进了现有文献的一些结论. 最后, 通过算例说明了该结果是可行的.
作 者：莫玉忠 谭伟冰 虞继敏 MO Yu-zhong TAN Wei-bing YU Ji-min 作者单位：莫玉忠,虞继敏,MO Yu-zhong,YU Ji-min(广西师范学院,数学与计算机科学系,广西,南宁,530023;柳州师范高等专科学校,数学与计算机科学系,广西,柳州,545004)
谭伟冰,TAN Wei-bing(广西师范学院,数学与计算机科学系,广西,南宁,530023)
混合线性模型 篇3
径流中长期预报一直以来是人们关注的重大课题, 在气候变化和人类活动共同的影响下, 径流的变化更加复杂化, 对其未来的精确描述十分困难。传统的径流预测方法主要有两大类:成因预测法和统计预测法[1]。成因预测法[2]多是基于研究大气环流、长期天气过程的演变规律和流域下垫面物理状况的确定性预测模型, 是径流预测研究的一个重要发展方向, 但离实用尚有较大差距。统计预测法[3]是从大量历史资料中应用数理统计的理论和方法, 寻找预报对象和预报因子之间的统计规律和关系或水文要素自身历史变化的统计规律, 建立预报模式进行预报, 常用的有时间序列法[4]、多元回归分析法[5]等。时间序列法具有原始资料搜集简便、预见期长等优点, 但它的缺点在于不能充分利用对预测量具有很大影响的气候信息和其他因素, 导致了预报的不确定性和数据的不稳定性。回归分析法虽然充分考虑了预测量对预测因子的依赖关系, 但它忽略了预测量自身的演变规律。本文将自回归和多元回归有机结合起来, 提出了一种非线性混合回归模型, 利用神经网络进行模型求解, 并应用于戴营站的年径流预报中。这种非线性混合回归模型既考虑了预测量自身演变的客观规律, 又利用了预测量对预测因子的依赖关系, 且充分考虑了水文过程的非线性关系, 模型具有一定的理论基础, 是一种值得研究的实用方法。
1 混合回归模型
1.1 多元回归模型
在影响中长期径流的众多相关因素中, 可以通过敏感性分析和弹性分析, 选择影响显著的因素, 建立多元线性回归预测模型 (Multi-Regression Model, MR) 。
设x1, x2, …, xp是p (p>1) 个线性无关的可控变量, y是表示径流量的因变量, 若它们之间的关系为:
则称为p元线性回归模型。式 (1) 中, b0, b1, …, bp, σ2都是与x1, x2, …, xp无关的未知参数, ε是随机误差。
自回归模型是在回归模型的基础上推广、导出的。
1.2 自回归模型
对于时间序列{yt}, 如果满足:
就称模型 (2) 为p阶自回归模型 (Auto Regression Model, AR) , 或者p阶自回归。式 (2) 中, a0, a1, a2, …, ap为自回归系数, ε1为模型相应的白噪声序列。
1.3 非线性混合回归模型
回归模型和自回归模型是统计模型中的两种最基本的模型, 尽管对许多随机现象, 用这两种模型都给予了很好的描述, 但对于复杂的水文现象, 单纯地使用某种模型去描述就表现出不足, 弥补单一模型不足的一个自然想法就是将两者结合起来就形成混合回归模型。文献[6]以年径流作为自回归因子, 降水、气温和用水作为回归因子, 建立了年径流预测的线性混合回归模型, 对黄河三门峡站年径流进行预测。对于理想的多输入、单输出系统, 设其输入为{x1}、{x2}、…、{x6}, 经过系统的作用后, 输出为{y}, 则线性混合回归模型为:
式中:yt为系统在t时刻的输出;a (0) , a (1) , …, a (s) 为系统的响应函数, 也称回归系数;q0, q1, …, qs为模型的阶数, 其中q0为自回归阶数, q1, …, qs为回归阶数;εt为模型的残差。
在全球气候变暖的背景下, 降水、气温等气象要素的变化引起径流的变化, 但是径流对降水、气温变化的响应并不是呈现简单的线性关系。刘昌明[7]通过大量坡面降雨试验从微观尺度上揭示了径流率随降雨强度的非线性变化现象。同时, 从宏观方面根据流域降雨径流观测资料, 采用流域水文模型的方法, 进行径流量对降雨和气温变化的响应的模拟分析, 同样证实年径流与年降水与年气温之间关系是非线性的, 影响径流变化的各因素之间是一种强耦合的非线性关系。
径流变化的影响因素与径流量之间是一个非线性映射问题, 因此, 为反应输入变量与输出变量之间的相关关系, 构造如下非线性混合回归模型:
式中:f (·) 称为非线性混合回归函数, 其中模型的阶数q0, q1, …, qs为待估计参数, q0为自回归阶数, q1, …, qs为回归阶数。
2 非线性混合回归模型的参数确定及预测效果分析
2.1 模型建立
径流形成过程是多种因素相互作用、相互联系的复杂的自然现象, 不同时间尺度下, 径流的变化受到的影响因子不一样。但总的来说, 在不考虑人类活动的情况下, 影响河流天然径流量变化的主要因素为气温、降水变化, 主要由于全球气候变化改变了与水循环密切相关的气候因子等;另外, 由于年径流的相依性, 年径流量自身也存在相关关系。因此, 流域水文系统可以看作是一个多 (因素) 输入、单 (径流) 输出系统。
在本文研究中, 输入因子为以往的年径流量{Qt-1, Qt-2, …, Qt-n}、降水量{Pt, Pt-1, …, Pt-n}和气温{Tt, Tt-1, …, Tt-n}, 输出因子为年径流量{Qt}, 由于径流变化的影响因素与径流量之间是一个非线性映射问题, 因此, 为反应输入变量与输出变量之间的相关关系, 构造如下非线性混合回归模型:
式中:f (·) 为非线性函数, q0, q1, q2是模型的阶数;εt为模型的残差。
上式 (5) 即为非线性混合回归系统水文模型的结构。
2.2 模型定阶
根据上述模型结构, 采用潮河流域戴营水文站 (1956-2000年) 的降水量、年气温和修正后的年径流来建立预测模型, 其中1956-1992年前37年的数据作为训练样本, 1993-2000年后8年作为预测检验样本。
以线性混合回归模型的阶数来确定非线性混合回归模型的阶数。线性混合回归模型的阶数确定和参数识别方法参见文献[7]。参数的估算和阶数的确定是相互联系、互为前提的, 可通过Matlab编程一次实现。计算得到自回归阶数为2, 回归阶数为2的线性混合回归预测模型如下:
因此, 影响未来径流的因子包括当年的降水、气温以及前期的降水、气温和径流量, 共有6个因子作为输入。
2.3 参数估计
根据Kolmogorov定理, 任何一个时间序列都可以看成是一个由非线性机制确定的输入输出系统的近似, 由于神经网络具有通过学习逼近任意非线性映射的能力, 故采用神经网络强大的非线性映射能力来实现式 (5) 的非线性函数f (·) 。
设p个样本集合{ (x (t) , y (t) ) |x∈Rm, y∈Rn, t=1, 2, …, p}, 式中R为映射空间。对于p个样本集合的离散空间序列, 神经网络可以完成从输入到输出的高度非线性映射, 即可以找到某种映射使得f:Rm→Rn。采用输入节点为m, 输出节点为n, 隐层节点为u的3层神经网络来实现f, 网络的输入与输出之间的关系如下:
式中:f为sigmoid函数;Φ为输入层至隐含层的激励函数;k=1, 2, …, n;t=1, 2, …, p;xi (t) 为t时刻网络的输入;
设网络总的误差小于ε, 即
根据kolmogorov定理[8], 确定非线性混合回归预测模型的BP网络隐含层节点数为13, 最终选择[1,6]网络结构进行年径流的预测。
BP网络的隐层神经元采用tansig函数作为激励函数, 输出层神经元采用logsig函数作为激励函数。采用自适应学习速率动量梯度下降反向传播traingdx函数对网络进行训练。规定的期望误差为0.01, 最多迭代次数为5000次, 学习效率为0.1。
2.4 数据处理
由于BP网络以值域为[0, 1]的S形函数作为转换函数, 因此在利用神经网络进行预测之前, 需将实际数据进行适当的变换, 规范到[0, 1]区间, 使数据主动适应神经网络, 以提高神经网络预测的准确性。文献[9]对于研究了线性变换和几种基本的非线性平滑变换, 指出对数据样本进行COS型变换的效果较好。本文采用COS型变换, 将数据样本变换到[0, 1]区间。COS型变换的基本形式为:
2.5 预测结果分析
为使训练后的网络具有较好的预测能力, 必须有足够的样本, 否则网络将无法归纳出样本集中的内在特征。但样本过多则会造成网络过度接近样本值, 从而丧失归纳和推理能力。经过多次实验, 最后选用1956-1992年共37年的数据作为训练样本, 以1993-2000年共8年的数据作为检验样本来检验模型的精度。按照上面的训练过程, 以matlab编程实现此算法, 得到1958-1992年的年径流量的拟合值如图1, 图2为戴营1958-1992年实测和模拟径流散点图。再通过网络训练得出的权值和阈值, 将后8年的检验样本输入此网络, 得到1993-2000年的各年径流预测值如图3, 图4为戴营1993-2000年实测和模拟径流散点图。
为了进一步说明非线性混合回归模型的优越性, 本文采用自回归模型 (Auto-regression Model, AR) 、多元回归模型 (Multi-Regression Model, MR) 和线性混合回归模型 (Linear Mix Regression Model, LMR) 对同样的算例进行对比计算, 4种方法的比较结果见表1。从表1所显示的结果可以看出, 基于BP网络的非线性混合回归模型的预报合格率为75%, 对径流预报效果较好, 误差变化范围缩小, 预报稳定且精度高, 充分说明了这种方法应用于径流中长期预报的有效性。
3 结 语
本文建立了以年径流量为自回归因子, 年降水和气温为回归因子的年径流预测的非线性混合回归模型, 利用BP网络进行模型求解, 并对戴营水文站的年径流进行预测, 拟合期的年径流计算值与实测值线性回归系数R2为0.825 4, 检验期的年径流计算值与实测值线性回归系数R2为0.747 3, 预测结果表明了非线性混合回归模型进行年径流预测的有效性。同时将非线性混合回归模型与自回归模型、多元回归模型和线性混合回归模型进行对比研究, 结果表明基于BP网络的非线性混合回归模型的预报合格率为75%, 对径流预报效果较好, 误差变化范围缩小, 预报稳定且精度高, 充分说明了这种方法应用于径流中长期预报的优越性。
参考文献
[1]Chow V T, Maidment D R, Mays L W.Applied Hydrology[M].New York:McGraw-Hill, 1988.
[2]Thompson S A.Simulating of climate change impacts on waterbalances in the central United States[J].Phys Geogr, 1992, 13 (1) :31-52.
[3]Freeze R A.A stochastic-conceptual analysis of rainfall runoffprocesses on a hillslope[J].Water Resources Research, 16 (2) :391-408.
[4]洪时中.非线性时间序列分析的最新进展及其在地球科学中的应川前景[J].地球科学进展, 1999, 14 (6) :559-565.
[5]王秀杰, 练继建, 费守明, 等.基于小波消噪的混沌多元回归日径流预测模型[J].系统仿真学报, 2007, 19 (15) :3 605-36-8.
[6]蒋晓辉, 刘昌明, 刘煜, 等.用混合回归模型预测黄河三门峡站年径流变化的初步研究[J].中国科学E辑, 2004, 34 (增刊) :95-102.
[7]刘昌明.黄河流域水循环演变若干问题的研究[J].水科学进展, 2004, (5) :608-614.
[8]张乃尧, 阎平凡.神经网络与模糊控制[M].北京:清华大学出版社, 2000.
混合线性模型 篇4
张璇 王嘉宇
2011-12-13 14:33:23 来源：《统计与决策》(武汉)2010年15期第4～8页
内容提要：文章运用Jackknife和Boostrap的方法，对参数估计的方差进行改进，构造了合适的参数估计的置信区间。通过样本组数和组内个体数的变化，利用数据模拟的方法进行研究，表明参数估计的可靠性很大程度上依赖于组数；对于固定效应参数，组数取30就可以得到可靠的估计值。对于σ和方差协方差成分T，组数分别取50和70才能得到可靠的估计。
关键词：分层线性模型 参数估计的覆盖率 Jackknife Boostrap 数据模拟
作者简介：张璇(1979-)，女，湖南湘潭人，中国人民大学统计学院博士研究生，讲师，研究方向：统计模型及其计算、计量经济学（北京100084）；王嘉宇，卡尔斯塔德大学国民经济与统计系，乌普萨拉大学信息科学与统计系（瑞典65188）。
1研究背景
很多社会研究都涉及分层数据结构，例如，经济学家探求在多个国家中经济政策是如何影响居民的消费行为，研究采集的观测数据不仅包括以国家为层次的经济指标，还包括以家庭为单位的信息，因此整个观测的数据结构是分层的。此时，同属一个层次的个体之间的相关性会大于来自不同层次的个体之间的相关性，整个观测样本就不再具有独立同分布性质，如果继续使用经典的线性回归模型，就会得到有偏的参数估计和错误的统计推断结果。近年来，随着分层线性模型统计理论的发展，一套完整的应用于分层结构数据的统计推断方法已经建立起来，并且能得到有效的参数估计。
分层线性模型(hierarchical linear models)的称谓最早由Lindley和Smith(1972)[1]提出。这个模型在不同的研究领域有不同的称呼，在社会学研究中，它经常被称为多层线性模型(multilevel linear model)；在生物统计研究中常用的名字是混合效应模型(mixed-effects models)和随机效应模型(random-effects models)；计量经济学文献称之为随机系数回归模型(random-coefficient regression models)等。
模型的一般形式为：
目前比较常用的分层线性模型的估计方法有极大似然法(ML)，约束极大似然法(REML)(Littell，1996)[2]和完全贝叶斯法(Full Bayes)。理论研究表明，大样本情况下ML和REML得到的参数估计都是一致最优估计量，但是在样本较小的情况下，REML在估计方差协方差成分确定性，因此REML比ML得到了更可靠的和T时，考虑了固定效应系数γ的不
和T的估计值。Full Bayes考虑
了所有辅助参数的不确定性，因此理论上得到的参数估计较ML和REML更可靠，但是此估计方法需要运用MCMC算法，很多情况下，并非能得到一个收敛的Markov链，当分层线性模型形式较复杂时（如待定参数向量增多、层数增加等），Full Bayes方法相当复杂。因此很多统计软件采用ML和REML估计分层线性模型，本研究采用REML方法估计。
由于ML和REML方法估计的前提假设都是样本量要足够大，因此样本量较小时，这两种估计都是有偏的，由此得到的参数置信区间和假设检验都是不可信的。因此，分层线性模型样本量问题的研究是一个重大课题，近十几年来，很多的学者都致力于这方面的研究。Bryk和Raudenbush(1992)[3]很早就指出，可以凭借OLS回归的经验法则：增加一个解释变量至少需要增加10个观测样本，联系到分层线性模型的估计，增加一个层2的结果变量（层1模型中待定的随机参数）至少需要增加10个观测样本，这个准则只是OLS回归经验法则的平移，并没有清楚回答模型要得到较好估计所需的最小样本量问题。Kim(1990)[4]在研究分层线性模型斜率参数估计的时候，发现当组数较少而组内的观测值相对很大时，固定效应和随机效应的估计会产生很大的偏差，但是Kim的研究在相同的条件下，只进行了50次模拟计算，因此结果不太可靠。Mok(1995)[5]的研究也得出了类似Kim的结论，并且他进一步指出参数估计的偏差和效率更大程度取决于样本中的组数大小。Busing(1993)[6]指出当组数超过300时，随机效应的估计才是无偏的。Kreft(1996)[7]运用模拟技术探讨了分层线性模型的势，建议30/30准则，即30个组，每组30个观测值可以得到较可靠的估计。Hox(1998)[8]在Keft(1996)之后更详尽阐述表明，50组，每组20个观测值可以使交互效应得到较可靠的假设检验；100组，每组10个观测值可以较可靠进行随机效应的假设检验。另外Bliese(1998)[9]明确指出组内相关系数(ICC)与样本量相互联系。Snijders和Bosker(1999)[10]阐明对于较大的ICC值，增加组数也不会得到一致的标准差的估计。Mass和Hox(2005)[11]针对不同的层1和层2样本量进行了模拟研究，结果表明层2的样本数小于50，会导致有偏的协方差成分估计。
纵观以上研究，其方法主要是运用数据模拟的方法，针对不同层
1、层2样本量的组合，估计分层线性模型，比较固定效应、随机效应和方差协方差成分的估计值。在计算机发展日新月异的今天，进行模拟计算简单易行，应用广泛。但是对于此类问题，以往研究存在许多问题：①大部分研究考虑的是参数估计的准确性，只有很少的一部分重视估计的标准差的准确性。②许多研究在讨论估计参数假设检验或构造置信区间时，想当然的运用了大样本的近似理论，即参数估计近似服从正态分布，但是实际上，在样本量较小时，分层线性模型的参数估计值，尤其是方差协方差成分的估计值并非服从正态分布(Raudenbush(1984))。
2研究方法
本研究仍采用数据模拟的方法，用R语言编译计算程序。另外，研究采用简单的两层模型，每层模型各有一个解释变量：
在数据模拟过程中，研究分别改变组数，组内个体数和组间相关系数ICC，在样本变化的条件下，构造参数估计的置信区间，计算置信区间覆盖真实参数的覆盖率来考查估计值的分布情况。①组内个体数取值为5，30，50，ICC取值为1，2，3，考察组数从5至100改变情况下，估计值的分布如何变化；②组数取值为30，50，100，ICC取值为1，2，3，考察组内个体数从5至100改变情况下，估计值的分布如何变化。
在构造置信区间时，以往研究直接将方差用近似标准差代替，构造正态分布的置信区间，而本研究用另外一种方法计算参数估计的标准差。
R程序lme4程序包中的lmer命令可以提供分层线性模型的计算，它提供固 定效应系数和方差协方差成分的估计值，以及固定效应系数估计的标准差，但是不提供方差协方差成分估计的标准差，lmer命令的编写者Doug Bates(2008)表示，第一，方差协方差成分估计的标准差的计算是非常复杂的，很多情况下，估计值已经在参数空间以外，即便是估计出来，也是无意义的，因此有些软件（如SAS）提供标准差也并不可靠。第二，当参数分布差不多是对称的情况下，提供参数估计的标准差才是有意义的，因为可以由此构造置信区间。
Harvey Goldstein(1999)提到，样本不是很大的情况下，如果随机扰动项不再服从正态分布，分层线性模型的固定效应和随机效应的估计仍然是一致的，但是其估计的标准差不能用来构建置信区间和进行显著性检验。一种替代的方法是运用密集计算中的Jackknife和Bootstrap计算估计的标准差，由Jackknife和Bootstrap理论可知估计的标准差可以利用正态分布构建置信区间和进行显著性检验。
3研究步骤
3.1模拟数据的生成(4)讨论覆盖真实值情况：如果以上的置信区间确实覆盖了真实值，则取值为1；否则取0。
(5)计算覆盖率：重复第一步至第四步1000次，统计覆盖次数，计算覆盖率。
(6)在不同样本量情况下，考察覆盖率随样本量如何变化。固定n(=30，50，100)，p(=1，2，3)，N从5到100变化；固定N(=5，30，50)，p(=1，2，3)，n从5到100变化。
3.2.2运用参数Bootstrap方法
和Jackknife方法相比，参数Bootstrap方法步骤(2)中有差异，它将(2)分成以下几步。
构造置信区间和覆盖率的计算和Jackknife的过程一样。
4研究结果
4.1固定效应和方差协方差成分的置信区间覆盖情况
由于不论是Jackknife还是Bootstrap方法，运算量都是非常大的，因此我们先考察在组成样本三个条件N，n，p取不同值的组合下置信区间的覆盖情况。
4.1.1Jackknife方法
表1的结果表明，当N=5时，计算的覆盖率都离95%较远，这时若n的取值较大，如n=100，覆盖率会相对提高，但是方差协方差成分中的覆盖率还是很低，这说明分层线性模型的估计很大程度上依赖于N的个数。如果仅仅考虑参数点估计，当N的取值很小，即使n取值很大，的估计值也不能收敛，因此存在某个整体较大的样本量（较大的n和较小的N），方差协方差成分的估计是相当不可靠的。另外从表1中，还可以看出组内相关系数(ICC)的改变对置信区间的覆盖率没有太大影响。
当N达到30，n达到30时，固定效应和的估计是比较可靠的，其覆盖率基本上达到93%，只有三个值小于93%；但是此时对于仍然偏低，只有个别能达到93%以上。
表1N，n，p不同取值的组合下，Jackknife和Bootstrap方法计算参数95%置信区间覆盖率(%)
注：“J”表示Jackknife法得到的覆盖率；“B”表示Bootstrap法得到的覆盖率。
4.1.2Bootstrap方法
而言，大部分的覆盖率
和Jackknife相比，Bootstrap方法得出的结论很相似，但是对于N较小的时候，可以清楚地看到ICC的值越大，其95%的置信区间的覆盖率越小；另外对于σ的估计和Jackknife不一样，当N=5时，只要提高组内的样本量n，例如当n=100时，其95%的置信区间的覆盖率的均值为94.5%；当n=30，N=5,95%的置信区间的覆盖率的均值也达到93%，这说明σ的估计依赖的是整个样本容量，若组数稍小时，可以通过增加组内的个体数，提高估计的可靠性。
从以上的分析可以看出，Jackknife和Bootstrap方法只是在样本组数较小时有差异，而且其差异主要表现在σ的估计上，这种差异来源于这两种方法计算时再抽样的过程不同，当然这种Parametric Bootstrap的方法在样本组数较小时优于Jackknife的方法，但当样本数增加时区别不大，而Parametric Bootstrap的计算量远远超过Jackknife的计算，因此本研究下面的分析，当n和N连续的变化时，我们考虑Jackknife的计算方法。
4.2N从5到100变化时参数估计置信区间的覆盖率
现在考虑n取30，50，100；ICC取1，2，3的不同组合，N从5到100变化。因为ICC的变化对覆盖率的影响不大，因此以下的覆盖率均为ICC三种取值下的平均覆盖率。另外，为了考察覆盖率的可接受性，我们利用随机模拟的方法建立相应样本量下的覆盖率的置信区间（利用随机数均值覆盖率的置信区间），构成了相应的置信带。
4.2.1固定效应
以为例，其他三个固定效应的变化情况相似。从图1清楚地看到，当N取值很小时（如N=5，10），固定效应的覆盖率是很低的；当N达到35时，三种情况下的覆盖率都进入了置信带。
4.2.2方差协方差
考虑σ和的变化参照
。①σ的估计情况；当N很小时（如N=5,10），标准差σ覆盖率远远低于置信带的下界；随着n的增加，有更多的覆盖率进入置信带，平均当N达到50时，三种情况下的覆盖率完全进入置信带。②计情况：的估的覆盖率在n的三种取值情况下，都比固定效应γ和标准差σ相对应的覆盖率低，因此N从5变化至100时，只有较少的部分进入置信带，对于n的三种取值，N达到70时，覆盖率几乎完全进入置信带。对于较小的N，随着n的增加，参数覆盖率更靠近置信带。因此对于τ的估计，可以通过n的增加得以改善。
4.3n从5到100变化时参数估计置信区间的覆盖率
现在考虑N取5，30，50；ICC取1，2，3的不同组合，N从5到100变化。因为ICC的变化对覆盖率的影响不大，因此以下的覆盖率均为ICC三种取值下的平均覆盖率值。另外，为了考察覆盖率的可接受性，我们利用随机模拟的方法建立相应样本量下的覆盖率的置信区间（利用随机数均值覆盖率的置信区间），构成了相应的置信带。
4.3.1固定效应
以为例，其他三个固定效应的变化情况相似。图4（见下页）清楚地表明，N为5时，不论n如何增加，参数置信区间的覆盖率都远离置信带；当N取30时，只有两点在置信带外；当N取50时，整条覆盖率连线全部进入了置信带。
图1n的不同取值下
图2n的不同取值下σ的95%的置信区间覆盖率
图3n的不同取值下
图4N取5，30，50时，4.3.2方差协方差 的95%的置信区间覆盖率 的95%的置信区间覆盖率 的95%的置信区间覆盖率
考虑σ和，的变化参照。①σ的估计情况。图5显示出，N为5时，不论n如何增加，参数σ置信区间的覆盖率都远离置信带；当N取30时，大部分的覆盖率都进入了置信带；当N取50时，除一点外，整条覆盖率连线全部进入了置信带。②的估计情况。图6表明的估计和固定效应参数与σ有很大差别，当N取50时，还有很多覆盖率在置信带以外，这表明对于方差协方差成分的估计N为50是不够的。结合图3，我们可以看到对于方差协方差成分的估计，N达到70或者更多才可靠。另外，图6还显示了一个有趣的现象，当N取30和50时，较小的n（比如n＜20）的覆盖率稳定地进入了置信带，而较大的n（例如：30＜n＜90）的覆盖率反而在置信带以外，这表明这时较小的n的参数估计比较大的n更可靠。
5研究结论
从以上的深入分析中，我们可以看到分层线性模型估计的可靠性与样本中的组数，组内个体样本数的取值有密切关系，总结如下：
(1)样本中的组数N。分层线性模型估计的可靠性很大程度上依赖于N的取值，随着N的增加，固定效应，方差协方差成分的参数估计的95%的置信区间的真实覆盖率能够进入置信带。当N较小时（如N=5），无论n如何增加，所有参数估计都不可靠，其覆盖率远远低于95%的置信带。对于固定效应参数γ的估计，N取30就可以得到可靠的参数估计。对于参数σ的估计，N取50，才能使其得到可靠的估计，方差协方差成分T的可靠估计需要N的取值达到70以上。
图5N取5，30，50时，σ的95%的置信区间覆盖率
图6N取5，30，50时，(2)组内个体样本数n。在N固定的情况下，随着n的增加，所有参数估计的可靠性都没有明显得到改善，即覆盖率连线没有明显的上升趋势。另外当N取值为30和70之间时，对于较小的n(n＜20)，比较大的n的参数估计更可靠。
参考文献：
[1]Lindley, D.V., Smith, A.F.M.Bayes Estimates for Linear Model[J].Journal of the Royal Statistical Society, Series B,1972,(34).[2]Littell, R., Milliken, G., Stroup, W., Wolfinger, R.SAS System for Mixed Models[M].Cary, NC: SAS Institute Incorporated, 1996.[3]Bryk, A.S., Raudenbush, S.W.Hierarchical Linear Models[M].Newbury Park, CA: Sage, 1992.[4]Kim, K.S.Multilevel Data Analysis: A Comparison of Analytical Alternatives[M].Los Angeles: University of California, 1990.的95%的置信区间覆盖率
混合线性模型 篇5
1 资料与方法
1.1 临床资料
选取2010年1月至2011年6月我院收治的癫痫患儿共180例，均按医嘱每日早晚连续服卡马西平剂量的时间均超过1个月，其血药浓度达稳态，根据治疗过程是否联合用药分为2组，实验组单用卡马西平，对照组采用卡马西平合并丙戊酸钠、癫痫宁片、癫痫平片、托砒酯、氯硝西泮等抗癫痫药中的1种或以上。两组患者的临床资料见表1。
按照临床的治疗需要，全部患儿均于早晨服药前空腹进行采集静脉血样及监测，并记录抽血前的末次服药剂量、服药时间、采样时间及患儿的体量、身高、生理及病理等状况。
1.2 方法
1.2.1 测定方法
采用荧光偏振免疫法测定全部患者的血清药物浓度，卡马西平稳态血药浓度结果见表1，各检测项目RSD%均<5%，检测均在各项目标注的线性范围内进行。数据采用NONMEM程序、Version VI、Level2.1双精度及SPSS10.0统计软件进行处理。
1.2.2 PPK模型的建立
药动学模型：
选择一级吸收与消除的一房室开放模型为药动学模型。公式：Cp=Dose×（Ka×K）×（e-Kt-e-Kat）/[CL×（Ka-K）]。其中Cp是模型预测的血药浓度；Dose是口服药物剂量；K是消除速率常数；Ka是吸收速率常数；CL是清除率。本研究收集稳态消除相的样品，对CL和表观分布容积进行估算。
固定效应模型：
基本模型不引入影响因素（即协变量），建立基本模型：CL=THETA(1）×EXP(ETA(1））;Vd=THETA(2）。CL是清除率；THE-TA(1）、THETA(2）分别是CL与Vd的群体典型值；ETA(1）是CL的个体间差异；Vd是表观分布容积。考察性别、年龄、体质量、服药天数、每日服药剂量及合并其他抗癫痫药物对群体药动学参数影响。分别将固定效应因素加入基本模型中，检验水平为0.05，如目标函数值改变>3.84，则保留该因素，反之则剔除。考察全部因素后，根据目标函数值下降的多少进行排序，依次从大至小将固定效应（协变量）引入基本模型，如OBJ的改变少于3.84则剔除该固定效应因素。
统计学模型。
个体间变异模型应用指数模型，公式为CLj=CL’×eηCL,j。CLj是第j个患者的相对CL;CL’是群体典型值，代表患者群体特性的相对CL；ηCL,j为当均数是0和方差是ω2时的个体间变异。个体自身变异模型（又称残差模型）应用加法模型，公式为Ci=C’ij+εij。Cij是第j个患者、第i个血药浓度的实测值；C’ij是相应的预测值；εij为当均数是0、方差是σ2时的残差变异。
1.3 模型的验证
采用内部有效性验证法中的自举法进行模型验证。自举法为目前模型验证流行使用的方法之一。若参数的频数分布接近正态时，以均数为点估计，应用正态原理进行估计可信区间；若参数的频数分布是偏态时，以中位数为点估计，以上下2.5%分位数作为95%置信区间。
2 结果
2.1 固定效应因素的初筛
根据观察在基础模型中分别引入各协变量后目标函数值下降的多少来判断对群体药动学参数的影响（表2）。
2.2 最终模型
根据目标函数值之差的大小，依次在基础模型中引入年龄、体质重及合并使用托砒酯，结果显示当引入托砒酯此项协变量后，目标函数值之差少于3.84，故剔除该因素。
2.3 终模型与基础模型比较
最终模型的群体预测值对个体预测值和观测值的特征图中的数据点集中趋势均强于基础模型，群体预测值与个体预测值及观测值的拟合程度均有很大的改善。最终模型更可靠、稳定。
2.4 最终模型验证
采用自举法进行模型验证，应用WFN提供的自举法检验抽样方法，次数B=200，并采用SPSS 11.0统计软件计算各参数的均值及95%置信区间。稳健率=成功数/抽样数×100%。结果显示稳健率高达96%，除体质重的θ值稍高于95%置信区间外，其他各参数值均位于95%置信区间范围中。
3 讨论
群体药动学（PPK）是由经典药动学结合统计学原理发展而来，使零散的血药浓度测定结果可用于群体参数值的估算，更加切合临床实际，为个体化给药的实施打下良好的基础[5]。非线性混合效应模型（NONMEM）[6]可对临床治疗时监测的稀疏数据进行数学分析，获得实际患者的药动学参数,优化给药方案。NONMEM法是迄今为止最被认可与采用的PPK参数测定法，在国外被广泛应用于个体化给药[7]。本研究通过建立卡马西平的群体药动学模型，揭示了对卡马西平在体内消除产生影响的主要因素及影响程度。本研究结果能为临床制定出更安全、有效、合理的卡马西平的给药方案。
参考文献
[1]中华人民共和国药典委员会编.临床用药须知是[S].2005:78.
[2]李金恒.临床治疗药物监测的方法和应用[M].北京:人民卫生出版社,2003:28.
[3]郑萍,胡敏燕,刘世霆,等.卡马西平吸收的多峰现象[J].中国临床药学杂志,1999,8(3):184-186.
[4]焦正,钟明康,施孝金,等.NONMEM法估算中国癫痫患者卡马西平的清除率[J].中国药理学通报,2003,19(3):314.
[5]张慧娟,李俊宏.高效液相色谱法测定儿童卡马西平药代动力学参数[J].临床药学,2004,13(3):64.
[6]焦正,钟明康,施孝全,等.NONMEM法考察中国癫痫患者丙戊酸钠和卡马西平的药动学相互作用[J].中国药学杂志,2004,39(2):130-132.
混合线性模型 篇6
过去几年里, 复杂网络得到越来越多的关注。很多系统都可以用复杂网络表示, 例如生物体的神经系统, 万维网, 代谢系统和电网等。这些网络的节点对应于系统的单元, 比如大脑中的神经元, 万维网中的超链接等, 而边则代表单元间的联系[1,2]。
在现实世界中, 网络的拓扑结构可能由于跳变或切换而发生变化。近年来, 带有切换拓扑结构的网络的同步性能得到了越来越多的关注[3,4,5,6,7,8,9,10]。切换经常使网络路径连接失败或者产生新路径。比如电网, 它是由很多发电机组成的一个网络系统, 只有当所有发电机都同步时, 电网才能正常工作。由多个局域电网组成的大型电网为了能正常供电, 必须在受到干扰的情况下也保持稳定。但是, 局部电网的不同步现象可能会影响整个电网, 甚至让整个电网瘫痪[4]。除了局部电网不同步之外, 新电网并入也可能破坏电网的同步性。为了解决这类问题, 一些新的切换拓扑网络的定理被提出。文献[5]证明了如果具有静态时间均值的拓扑结构的振荡器网络能够同步, 则当时间均值足够快的时候, 具有时变拓扑结构的网络同步。文献[6]用快切换方法研究了有向动态切换网络的自适应全局同步问题。由于控制系统的原件不能承受如此高的频率, 因此快切换在实际应用中无法推广。
慢切换方法包括共同李雅普诺夫函数法, 多李雅普诺夫函数法, 平均驻留时间法等。共同李雅普诺夫函数可用于证明复杂切换网络的同步性[7], 但是只有当所有子网络都对共同李雅普诺夫函数有解时, 这种方法才适用。当任何一种单独的拓扑结构都不能使切换网络同步时, 可以采用多李雅普诺夫函数处理这类问题。文献[11]通过多李雅普诺夫函数设计了一系列切换定律使得复杂动态网络同步。驻留时间法是解决切换网络同步问题的另一种方法。文献[9]使用ei (t) =xi (t) -si (t) 证明同步性问题, 其中在该证明中需要用到f (s (t) ) 在s (t) 处的雅克比矩阵, 因此f (xi (t) ) 在s (t) 处需可微, 这在实际应用中有很大局限性。文献[9]描述的模型是无时延系统, 时延只存在网络拓扑中, 但在实际情况中, 很多网络系统由于信号传输的延迟导致系统时延。文献[12]对具有网络时延特性的网络控制系统 (NCS) 进行了稳定性分析。本文采用了具有子系统时延和节点间耦合时延的双时延系统, 研究了线性混合耦合切换网络的同步性问题并且解决了f (xi (t) ) 不可微的限制问题。基于这个网络模型, 利用李雅普诺夫函数、平均驻留时间和线性矩阵不等式 (LMIs) , 推导出了一个时延相关的全局指数同步的充分条件, 并通过一个实例验证了结论的正确性。
1 系统描述
为了表示现实生活中的一些拓扑结构会变换的网络, 我们考虑一个带有双时变时延的线性混合耦合切换网络:
其中i=1, 2, …, N, C=diag (c1, c2, …, cn) ∈Rn×n是对角阵且ci>0, i=1, 2, …, n。A= (aij) n×n和B= (bij) n×n分别是权重和时延权重矩阵, I (t) = (I1 (t) , I2 (t) , …, IN (t) ) ∈Rn×n是外部输入。G= (Gij) N×N是网络拓扑结构的耦合构造矩阵, 其中Gij定义如下:如果节点i到节点j有连接, 则Gij=Gji=1, 否则Gij=Gji=0 (j≠i) , G的对角元素定义如下:
D= (dij) ∈Rn×n和Dτ= (dτij) ∈Rn×n是内耦合矩阵, f (xi (t) ) (i=1, 2, …, n) 是R上的单调不减函数。σ (t) :[0, ∞]→χ={1, 2, …, m}为切换信号, 是一个右连续的分段常函数, 0≤τ (t) ≤r。以下给出一些有用的定义和引理。
定义1定义同步集合S={ (x1T, x2T, …, xNT) T∈Rn×N:xi=xj, i, j∈N}。
定义2[3]对于任意T2>T1≥0, 记Nσ (t) (T1, T2) 为 (T1, T2) 上切换信号σ (t) 的个数。若当Ta>0, N0>0时有Nσ (t) (T1, T2) ≤N0+ (T2-T1) /Ta, 则称Ta为平均驻留时间。
假设1 fi (xi (t) ) (i=1, 2, …, n) 为李普希茨连续, 即存在常数Fi>0满足:
假设2τ (t) 是t的有界可微函数, 满足0≤τ·≤h<1, 其中h是正实数。
引理1[8]克罗内克积有如下性质:
引理2[10]若N×N阶矩阵G满足 (2) , 则存在 (N-1) × (N-1) 阶矩阵H满足MG=HM, 其中H=MGJ, M和J定义如下:
引理3[11]与以下任意条件等价:
其中
引理4[13]假设向量方程ω:[0, r]∈R满足下列积分。对任意对称阵W∈Rm×m和标量r>0, 有:
2 主要结果
本节基于李雅普诺夫函数, 平均驻留时间和LMIs, 推导出一个关于复杂网络式 (1) 的全局指数同步的充分条件。
事实上, 切换信号σ (t) 和一系列拥有相同输入和不同的内耦合及外耦合矩阵Dk和Gk的子网络式 (6) 构成了切换网络式 (1) 。
为了表述简单, 记AфB为矩阵A和B的克罗内克积。
其中IN是N阶单位阵。则复杂网络式 (1) 可重写成:
定理1对于一个给定的常数α>0, 若满足假设1和2, 且fi (xi (t) ) 单调不减, 对于任意满足平均驻留时间的切换信号, 其中μ>1满足其中l, k∈χ, 为 (N-1) n× (N-1) n阶正定矩阵, 有 (N-1) n× (N-1) n阶正定对角阵使得, 其中:
其中M和J在式 (4) 中定义。则复杂网络式 (7) 能达到全局指数同步。
证明考虑下列李雅普诺夫函数:
其中:
在每个区间[tr, tr+1) 上的切换信号σ (t) 是常值。为了简便, 当t∈[tr, tr+1) 时记σ (t) =k, 此时对V1k微分得:
根据引理1可得:
由此可得:
由假设2知, 且τ (t) ≤r, 对V2k微分得:
同理得到V3k的微分:
对V4k微分得:
由 (5) 可得:
设根据式 (9) -式 (13) 得到:
根据式 (3) , 可得:
根据式 (14) -式 (18) 得到:
其中:
根据引理2得其中定义在定理1中。因此, 若则
在区间[tr, t) 上对积分得:
若存在一个常数μ≥1满足使得在切换时刻, 即t=tr, 有:
根据式 (20) 、式 (21) 和q=Nσ (t) (t0, t) ≤N0+ (t-t0) /Ta得到:
根据V (t) 的定义可得:
其中‖·‖是欧几里得范数, 其中分别是的最大特征根, 是最小特征根。
结合式 (22) 、式 (23) 可得:
因此当时复杂切换网络式 (1) 达到全局指数同步。证明完毕。
3 数值仿真
为了验证定理的正确性, 文章给出了下面的例子。为简单起见, 假设有两个子网络, 即χ={1, 2}, 每个子网络有10个节点, 即N=10。假设τ (t) =0.02[1+sin (40t) ], 则为了满足假设2, 取r=0.04, h=0.8。例子:考虑一个二维时延系统作为一个节点, 描述如下:
显然F=I2符合假设1。设外耦合阵分别为:
内耦合矩阵和时延内耦合矩阵分别为:
根据定理1, 通过解LMI不等式Ω~<0, 得到使这两个子网络达到全局指数同步的最佳衰减指数分别为α1=69.365, α2=5.843。取α=min{α1, α2}=5.843。解定理1中的LMI不等式得到μ=81.29, 根据定理1可知当切换信号的平均驻留时间Ta>Ta*时, 复杂切换网络式 (1) 达到全局指数同步。图1表示Ta=2 s时, 切换网络式 (1) 的误差曲线。图2表示Ta=2 s的切换信号。这两幅图说明当切换信号的驻留时间Ta=2 s>Ta*=0.75 s时, 切换网络可达到全局指数同步。图3表示Ta=0.5 s时, 切换网络式 (1) 的误差曲线, 图4表示Ta=0.5 s的切换信号。图3和图4表明当切换信号的驻留时间Ta=0.5 s
4 结语
混合线性模型 篇7
近几十年来磁悬浮技术得到了飞速发展, 而悬浮控制技术是磁悬浮技术中的核心和关键。悬浮控制系统性能的好坏, 将直接影响到磁浮系统的稳定性、安全性。为了达到悬浮的稳定性, 必须采用控制装置, 可以利用悬浮气隙的反馈对电流进行控制。磁浮系统的自适应控制、模糊控制、滑模变结构控制、神经网络控制以及鲁棒控制等方法都被深入地研究过。一般的控制方式是将非线性磁悬浮模型在平衡点处进行线性化, 然后对该线性化模型进行反馈控制。
一般线性控制基本上是利用气隙反馈量得到气隙变量, 由 PID调节后输出驱动电磁铁来达到调节气隙的目的。本文在一般线性控制方法的基础上进行了改进。利用简化的混合单电磁铁的动态模型得出数学模型, 而后在一般控制方法基础上得出以电流作为控制输入变量的控制模型。对一般线性控制方法和改进后的控制方法进行仿真比较分析, 得出了改进后的控制方法优于一般控制方法的结论。
2 混合单电磁铁悬浮模型
为了便于分析, 将混合悬浮模型简化为图1所示。图1中, 悬浮气隙的长度为z, 等效磁体的横截面积为A, 磁极铁芯与铁磁轨道中磁路的总长度为lm, 铁心的导磁率为μFe, 真空导磁率为μ0, 永磁体的厚度为l, 设永磁材料的剩磁为Br, 剩余矫顽力为Hc。永磁体的相对磁导率为μr=Br/ (HCμ0) 。
现假设:1) 忽略漏磁通;2) 磁势均匀降落在气隙上;3) 电磁铁仅有垂直方向上的移动, 其他方向受限无运动。则磁路方程为
式中:B为磁路中的磁感应密度。
由式 (1) 可得
令:
式中:I0为将永磁体的作用折算为同匝数为N的电磁铁时的等效励磁电流;z0为磁路中的永磁体以及铁磁材料折算过来的等效气隙的值。
则式 (2) 变为
式中:ie为将混合悬浮磁体等效为纯电磁铁时的励磁电流;ze为将整个磁体等效为纯电磁铁时的等效气隙。
磁体与轨道之间的吸引力大小为
根据牛顿第二定律可得系统的力学方程为
式中:m为悬浮磁极的质量;g为重力加速度 (9.8 m/s2) ;fd为扰动力。
假定在气隙为z0时稳定悬浮所需的励磁电流大小为i0, 则有如下关系成立:
设加在磁极励磁绕组两端的电压为u, 励磁绕组的总电阻为R, 则有:
式中:Ψ为任意时刻磁链;L (ze, ie) 为绕组中的电感。
将式 (5) 和式 (7) 在工作点 (i0, z0) 处线性化并结合式 (6) 的边界条件可得:
式 (8) 中的变量均为相对于平衡点处的增量, 常数L0, Ki, Kz分别为气隙电感、电流常数和气隙常数, 其表达式如下:
3 模型仿真与分析
模型仿真参数见表1。
3.1PID控制时的系统闭环结构图
引入状态变量
根据式 (9) 可得仿真模型框图如图2所示。
3.2仿真波形
气隙给定2 mm, 悬浮体质量分别为m=30 kg, 50 kg, 100 kg时的气隙响应仿真图如图3所示。
从图3a的仿真波形可以看出, 气隙响应开始处出现超调, 约为0.8 mm;图3b中超调量约为1.1 mm;图3c中超调量约为3.3 mm。由图3可以看出, 当气隙给定时, PID系统可以在极短的时间内使系统恢复稳定, 但都有一定的超调量, 且超调量大小随着悬浮体质量增加而增大, 由此造成电流的迅速增大。
4 模型改进
鉴于以上分析可知, 由式 (9) 所得的状态模型主要是悬浮启动超调量较大, 原因是开始时的电流误差太大, 导致PID调节时的开始电流增加太大所致。如果用一种梯进式电流增加方式给定, 不但使电流增加梯度很小、近于平滑, 而且也会使得悬浮气隙增加量较小, 达到悬浮稳定;待到悬浮稳定后再切入系统PID控制调节方式, 这样就不会出现开始悬浮时的超调过大现象。悬浮启动的控制框图如图4所示。
图4中, z0为初始气隙, i0和i1为梯进电流值, 初始悬浮气隙改变之前用梯进电流i0, 改变后用梯进电流i1, 且i1<i0, 以保证悬浮时平稳。z1为稳定气隙值。改进后的悬浮启动仿真框图如图5所示。
现以参数z0=5 mm, z1=3 mm, m=30, 50和100 kg分别进行仿真, 得出电流梯进和气隙变化曲线图如图6所示。
图6为m=100, 50, 30 kg时的电流与气隙变化曲线, 从中可看出:气隙从5 mm到3 mm之间变化缓慢;前段5 mm保持不变是因为电磁吸力小于悬浮体重量, 导致悬浮体无法悬浮;电流近似直线变化, 在前段5 mm保持不变时变化率较大, 后端5 mm到3 mm悬浮时较小;开始时悬浮力小于悬浮体重量, 导致悬浮体无法悬浮;电流近似直线变化, 悬浮时间随着悬浮体重量增加而延长。
在悬浮启动得到平稳基础上, 切换到PID系统调节控制, 这样可实现悬浮控制的平稳性, 图7为悬浮稳定后切换到系统PID调节控制时的气隙响应曲线, 由图7中可看出, 切换时出现了很小的超调, 相比PID单独控制时的超调量要小得多。
5 结论
将混合悬浮磁体简化为单混合电磁铁进行建模, 由此导出其状态空间方程;对一般的悬浮PID控制方式特性进行了仿真, 主要缺点为悬浮启动时的超调量较大。在此基础上提出了一种梯进式电流增量悬浮启动方法, 并对其进行了仿真分析, 验证了该方法能够提高悬浮体启动时的稳定性, 为更好的研究悬浮控制系统提供了有价值的参考资料。
参考文献
[1]Chan Mei-lang, Tao C W, Lee Tsu-tian.Sliding ModeController for Linear Systems with Mismatched Time-var-ying Uncertainties[J].Journal of the Franklin Institute, 2000, 337 (3) :105-115.
[2]Shafai B, Beale S, Larocca P, et al.Magnetic Bearing Con-trol Systems and Adaptive Forced Balancing[J].ControlSystems Magazine, IEEE, 1994, 14 (12) :4-13.
[3]Jie Xiao, Kulakowski B T.Sliding Mode Control of ActiveSuspension for Transit Buses Based on a Novel Air-springModel[C]∥American Control Conference, 2003, 5:3768-3773.
[4]Barr A J, Ray J I.Control of an Active Suspension UsingFuzzy Logic[C]∥Fuzzy Systems, 1996.Proceedings of theFifth IEEE International Conference, 1996, 1:42-48.
[5]Chiou Kuo-ching, Huang Shiuh-jer.An Adaptive Fuzzy Con-troller for 1/2 Vehicle Active Suspension Systems[C]∥Sys-tems, Man and Cybernetics.IEEE International Confer-ence, 2006, 2:1010-1015.
[6]Jin Yao, YuDejie, Song Xiaolin.An Integrated-error-basedAdaptive Neuron Control and Its Application to VehicleSuspension Systems[C]∥Control and Automation, 2007.ICCA 2007.IEEE International Conference, 2007, 1:564-569.
[7]Fu L J, Cao J G, Liao C R, et al.Study on Neural Net-works Control Algorithms for Automotive Adaptive Sus-pension Systems[C]∥Neural Networks and Brain, 2005.ICNN&B 2005.International Conference, 2005, 3:1795-1799.
[8]AL-holou N, Dae Sung Joo.Fuzzy Logic Based SuspensionSystem[C]∥System Theory, 1994.Proceedings of the26th Southeastern Symposium, 1994:172-176.
[9]谢云德, 常文森.电磁型 (EMS) 磁悬浮列车系统铅垂方向的建模与仿真[J].铁道学报, 1996, 18 (4) :47-54.
[10]胡基士.EMS型磁浮列车悬浮力分析[J].西南交通大学学报, 2001, 36 (1) :44-47.
混合线性模型 篇8
工业过程控制中,不同程度地存在时间滞后问题,较大的时滞会影响系统的稳定性,往往会使系统出现较大的超调量和较长的调节时间,甚至会出现震荡、发散,系统的控制性能明显变差。专家学者为解决这类问题作了有益地探索［1 － 5］,产生了许多优秀的控制算法( 例如Simth预估控制和智能控制等) ,但在这些方法中大多数存在着对被控对象精确数学模型要求高等缺点。基于此,本文将线性预测控制和经典PID控制有机结合形成一种复合控制方法,在一定程度上克服了上述方法的缺点。
本文提出用线性预测模型解决上述问题的方法是利用系统过去的几个输出值预测系统未来的输出值,将此预测值与期望设定值进行比较得出的偏差作为PID控制的输入,依照PID控制律来设定控制器的输出,从而使被延迟了的被控量超前反映到控制器,使控制器提前动作,实现事先调节,从而减少超调量和调节时间,消除时滞对系统控制品质的影响。该预测控制的优点是不依赖被控对象的数学模型,运算量小且响应速度快,适合复杂工业过程控制。另外,为了得到PID的理想控制效果,本文采用混合粒子群算法对预测模型系数和系统参数进行了在线优化和在线调整。实例仿真结果表明了本文提出的复合控制方法和优化方法远好于传统控制方法,大幅度减少了系统的超调量和调节时间,提高了系统的静动态性能,增强了系统的稳定性,全面提升的系统的控制品质。
1 线性预测模型
线性预测就是在假设已经知道y[t - 1],y[t - 2],…y[t - p]的基础上,最佳地预测当前尚未观察到的信号样本y[t],即把当前信号样本的估计值表示为过去p个信号样本的线性组合。
最佳预测就是要通过选取适当的预测系数使得式( 1) 中定义的预测误差e[t]的功率最小:
这里E[·]表示数学均值。
p的最优值可由AIC准则［6］确定。
未来d步的预报值为:
1) d = 1
2) 1 < d≤p时:
3) d > p时:
其中d=1,2,…为预报步数。当阶数p确定,θ(t)已知时,便可利用(3)~(5)式对系统未来输出y[t+d]进行预报,预报值为
这里指出,如果阶数p过大,会使运算量和储存量较大,影响实时性。
2 基于线性预测模型的PID控制系统
工业过程控制中的时滞现象表现为,给控制对象加一个控制量后,要经过一个较长的滞后时间d才能看到控制效果,致使系统产生明显的超调,使得系统的稳定性变差,调节时间延长,控制难度加大。利用线性预测模型的超前预测功能,对系统未来的行为进行预测,提前预测出系统的变化趋势,并与PID相结合,可在一定程度上克服时滞对系统控制造成的不利影响,使控制品质得以改善。控制系统结构如图1 所示。
该系统是用时刻t,t - 1,t - 2,…,t - p的输出y[t],y[t - 1],…,y[t - p]来预测未来d步的输出将此预测值作为反馈信号与期望设定值进行比较得出偏差,作为PID控制的输入,依照PID控制律来设定控制器的输出,从而使被延迟了的被控量超前反映到控制器,消除时滞对系统控制品质的影响。
3 基于混合粒子群算法的模型系数优化
PID控制器的一般形式:
PID的离散形式[7]:
式中k为采样序号,T为采样时间。
PID控制器的性能取决于它的三个参数( 比例Kp、积分Ki和微分Kd) 的整定是否合理,因此,优化PID控制器参数具有重要意义。但是PID参数的整定一般需要经验丰富的工程技术人员来完成,这种方法不仅费时,而且不能保证获得最佳性能。混合粒子群算法( HPSO) 是一种能有效避免早熟且具有很强的全局寻优能力的优化算法,本文将使用HPSO对PID控制器的参数和预测模型的系数进行优化设计,提高了系统的收敛速度。
3. 1 混合粒子群算法
标准粒子群算法通过追随个体极值和群体极值完成极值寻优,虽然操作简单,且能够快速收敛,但随着迭代次数的不断增加,在种群收敛集中的同时,各粒子也越来越相似,可能在局部最优解周边无法跳出。混合粒子群算法［8］摒弃了传统粒子群算法中的通过跟踪极值来更新粒子位置的方法,而是引入了遗传算法中的交叉和变异操作,通过粒子同个体极值和群体极值的交叉以及粒子自身变异的方式来搜索最优解。基本算法步骤如下:
1随机初始化种群中各微粒的位置和速度;
2评价每个粒子微粒的适应度,将当前各微粒的位置和适应值存储在各微粒的Pbest中,将所有Pbest中适应值最优个体的位置和适应值存储于Gbest中;
3更新每个微粒的速度和位置;
式中c1和c2为学习因子,w为常惯性权重。r1和r2为0 到1 之间均匀分布的随机数。
4 对每个微粒,将其适应值与其经历过的最好位置作比较,如果较好,则将其作为当前的最好位置;
5比较当前所有Pbest和Gbest的值,更新Gbest;
6将整个粒子群按适应值排序,用群体中最好的一半的粒子速度和位置替换最差的一半的位置和速度,保持Pbest和Gbest不变;
7若满足停止条件( 通常为预设的运算精度或迭代次数) ,搜索停止,输出结果,否则返回3继续搜索。
3. 2 实例分析
设被控对象为一延迟对象［9］:
采样时间为20 s,延迟时间为4 个采样周期,即80 s,被控对象离散化为:
PID控制参数kp= 0. 80,kd= 10,ki= 0. 002[10],取步长d =1,模型阶数p = 5 进行仿真实验。
5 阶线性预测模型为:
式中ai∈[- 2,2]( i = 1,2,3,4,5)
输入信号为一阶跃信号r( t) = 1。在混合粒子群优化算法中设惯性因子w = 0. 7,加速常数c1= 2,c2= 2,维数为5( 5 个待优化的模型参数a1,a2,a3,a4,a5) ,粒子群规模为30。为了提高控制系统的整体性能,将目标函数取为［11］:
式中( 读者可将目标函数取为式( 2) 进行实验) 。粒子算法的最大迭代次数为100,速度范围为[- 1,1],5个待优化参数范围均为[- 2,2]。利用MATLAB编程,运行得到参数的优化结果为:
仿真结果见图2。
实验表明,随着模型阶数的增加系统的超调量越来越少,调节时间越来越短,稳定性越来越高。但随着模型阶数的增加,运算量增大,需要更多的储存空间,时效性会变差。本文认为本算例取p=5较理想。
作为有效性的对比,下面采用经典PID算法对算例( 9) 进行仿真( 系统参数不变,即kp= 0. 80,kd= 10,ki= 0. 002) ,仿真结果见图3。
可以看出,基于混合粒子群算法的线性预测PID控制具有良好地静动态特性,能够大幅度地缩短系统的调节时间,使系统快速收敛; 另外,还强有力地抑制了系统的超调量,实现了参数的在线调整,提高了系统的控制品质,效果远好于传统PID控制,为解决大时滞问题给出了一种新思路和新方法。
另外、可以采用混合粒子群算法对PID系统参数kp,kd,ki及线性预测模型系数同时优化,相关参数设定同上。p = 5 时得到的优化参数值为:
仿真结果见图4。
从图4可以看出控制效果并没有得到有效改善。所以,要提高具有延时问题的控制效果及控制品质,除了对系统参数优化整定外,更应建立在对预测控制方法的改进和优化上,只有这样才能从根本上解决大时滞问题。
4 结束语
混合线性模型 篇9
考虑如下伪抛物方程的初边值问题
式(1)中为Hamilton算子，ΩR2为有界区域。p(x,y)，q(x,y)，F(x,y,t，·)，
f(x,y)为已知函数，且满足以下条件:
(2)，且满足连接性条件f(x,y)=0,(x,y)∈Ω;
(3)F为连续函数，且关于各变量存在连续的一阶偏导数。
涂慧等人在文献[1]中提出了一类非线性伪抛物方程的预测校正格式，曹志等也在文献[2]中提出了一维非线性伪抛物方程的混合体积元方法。本文研究二维问题(1)的混合体积元方法。
1 混合体积元格式
。则式(1)改写为
引入函数空间
从而问题式(2)的弱形式为:
求{u，σ}:(0,T]→W×V。使得
作区域Ω的拟一致三角剖分Th={K}，其中h=Km∈aΤxh{diam K}，在剖分Th的基础上作区域Ω的重心对偶剖分Th*={K*}，对偶网格由内部四边形与边界三角形组成，部分单元如图1所示。
图1中pi(i=1,2,3,4,5)为三角剖分的部分节点，对于内部节点p3，取B1,B2分别为ΔA1A2A3，ΔA1A3A5的重心，连接A1,B2,A3，同时连接A1,B1,A3，则四边形A1B1A3B2为相应于节点p3的对偶单元。
对于边界节点p5，取ΔA3A4A5的重心B3，连接A4,B3,A5则ΔA4B3A5为相应于节点p5的对偶剖分单元。
定义VhV为Th上的最低价的RaviartThomas空间，Wh为Th上的分片常向量空间。引入迁移算子γh，对vh∈Vh有
其中χQ为集合Q上的特征函数，N为原始剖分Th的边的总数目。
令Yh=γhvh，则上述问题的混合体积元格式为:求(σh,uh):(0,T]→Vh×Wh，使得
式(4)中uh(0)取为f(x,y)的混合体积元投影(见第3部分定义)。
2 一些引理
我们给几个重要的引理。
引理1[3]γh为自共轭算子，即
特别地，当φh为分片常向量时，
存在与h无关的正常数C和C0满足
引理2[3]存在与h无关的常数C满足
3 混合体积元椭圆投影
为了数值分析，引入如下的混合体积元投影。设满足
类似于文献[3]中，我们可以证明一下引理。
引理3当h充分小时，存在不依赖于h,t的常数C使得
4 连续时间的误差估计
将式(4)减去式(3)，并考虑到引入的混合体积椭圆投影，可得方程
令
我们可以将式(6)变换为
定理1设(u，σ)，(uh，σh)分别是式(3)和式(4)的解，并且满足所需要的正则性假设，则有如下误差估计
证明在式(7)中，令vh=ξ，wh=e+et，两式相加得，
下面我们分别估计式(8)的左右两边。对于左端有
下面分别处理式(8)右端各项。
从而得
当ε充分小时,
由Gronwall不等式有
则
由三角不等式可知定理结论成立。
参考文献
[1]涂慧,刘超,江成顺.二阶半线性为抛物方程差分格式与数值模拟.工程数学学报2,006;23(1):187—190
[2]曹志.两类拟线性发展方程的数值分析.硕士论文.济南:山东师范大学,2010
【混合线性模型】推荐阅读：
线性混合随机效应模型05-14
非线性混合效应模型法05-13
多元线性模型06-01
多元线性回归预测模型06-25
非线性函数模型07-29
混合模型05-12
混合多模型08-13
非线性组合预测模型08-21
混合DEA模型06-30
混合模型预测控制08-29