数字精度PID控制器改造的数值程序转换在数值精度、程序验证中的应用

可在www.sciencedirect.com在线获取理论计算机科学电子笔记317（2015）47-54www.elsevier.com/locate/entcs数字精度PID控制器的改造N. Damouche 1，CIM. 马特尔2号，阿拉巴马州。Chapoutot3，拉克鲁斯佩皮尼昂大学ViaDomitia，DALI，佩皮尼昂，法国法国蒙特佩利耶大学法国国家科学研究中心，LIRMM，UMR5506，蒙特佩利耶，法国巴黎科技大学，法国摘要执行浮点计算的数值程序对公式的编写方式非常敏感。为了提高表达式的准确性，已经提出了几种关于表达式转换的技术，现在我们的目标是通过自动转换更大的片段更进一步包含多个赋值和控制结构的代码。本文提出了一个案例研究在这个方向。我们考虑一个PID控制器，我们转换其代码，以提高其精度。当我们比较代码的不同版本（数学上是等价）表明这些变换对计算的准确性有重大影响关键词：数值精度，基于语义的程序转换，浮点运算，数值程序验证。1引言执行浮点计算的数值程序对公式的编写方式非常敏感。事实上，算术表达式中的小的语法变化并没有改变它们的数学意义，这可能会导致它们的求值结果发生重大变化。这种对表达式书写方式的敏感性是由于浮点运算的特殊性，例如，加法不是结合的，乘法不是可逆的[1，11，12]。此外，很难猜测公式的哪种写法在使用浮点数进行评估时具有最佳精度。这些年来，抽象解释技术[3，5]已经发展到推断计算结果的舍入误差的安全近似[2，6，7]。1电子邮件：nasrine. univ-perp.fr2电子邮件：univ-perp.fr3电邮地址：alexandre. ensta-paristech.fr†这项工作得到了ANR项目ANR-12-INSE-0007“CAFEIN”的支持http://dx.doi.org/10.1016/j.entcs.2015.10.0061571-0661/© 2015作者。出版社：Elsevier B.V.这是一篇基于CC BY-NC-ND许可证的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）。48N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）47我们的工作涉及到自动转换的计算，以提高其数值精度。关于表达式的转换已经得到了几个结果[9，10]，现在我们的目标是通过自动转换包含多个赋值和控制结构的较大代码片段更进一步。本文提出了一个案例研究在这个方向。我们考虑一个PID控制器，我们转换其代码，以提高其精度。更准确地说，我们采取一个初始PID代码，我们应用于它的几个处理，以产生其他PID程序，这是数学上等价的初始一个，但更准确的计算。第一个转换只重写赋值，而在第二个转换中，循环被展开。虽然这些转换是手工完成的，但它们被系统地应用，我们的目标是在未来的工作中自动化。当我们比较三种代码（数学上等价）的执行时获得的实验数据表明，这些重写对计算的准确性有重大影响本文的其余部分组织如下。第2节介绍了原始控制器PID1。在第3节和第4节中进行了转换，得到PID2和PID3。实验结果在第5节和第6节总结。2 PID控制器在本节中，我们将简要描述清单1中的原始PID程序。这种算法用于嵌入式和关键系统中，以将测量值m保持在称为设定点c的特定值。误差是设定点和测量值之间的差，控制器根据误差的积分i和导数d以及比例误差项p.计算这些项的加权和。权值用于提高程序的反应性、鲁棒性和速度这三个术语是：i) 比例项p：误差e乘以比例因子kp，p = kp×e。ii) 积分项i：误差e被积分并乘以积分因子ki，i= i+（ki× e × dt）。iii) 导数项d：误差e相对于时间被微分，然后乘以导数因子kd。让eold表示上一次迭代的误差值，我们有：1d= kd ×（e-eold）× dt。在实践中，存在许多其他方法来计算项d和i。在我们的实现中，它们是通过欧拉方法和矩形方法计算的N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）4749ni=0kp+kd×dt+i0+ki×dt×e i− k d× e n−1× dt。（i=1清单1：PID1的源代码。kp= 9.4514;ki = 0.69006;kd = 2.8454;invdt =5.0;dt=0.2;m= 8.0;c = 5.0;eold= 0.0;而真做e=c- m;p=kp*e;i= i+ki*dt*e;d= kd * invdt *（ e-eold）; r=p+i+d;eold= e;/* 更新内存*/m=m+0.01*r;/*计算测量：计划t*/分别下一节中描述的转换与这些特定算法无关。3、如何获得更精确的PID？在本节中，我们将详细介绍将原始PID（名为PID1）转换为新的等效程序（名为PID2）所需的不同步骤。 其主要思想包括开发和简化PID1的表达式，并将其内联到循环中，以提取常量表达式并减少操作次数。在循环的第n次迭代中，我们有：pn=kp×en，in=in−1+ki×en×dt，1d n= k d×（e n−e n−1）×dt。如果我们将pn，in和dn的表达式内联到结果表达式的公式中，在提取公因子后，我们得到：.1ΣΣ1然后我们注意到方程（1）中存在一些常数子表达式。所以，我们在进入循环之前计算一次我们有：1 1c1= kp+ kd× dt，c2 = ki× dt，c3 = kd×dt。接下来，我们在变量s中记录和s=n−1ei，它将不同的误差从e0到en-1。 最后，我们有：r n= R + i n，其中 R = c1× e n− c3× e n−1。我们的PID2算法如清单2所示。rn=en×50N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）47i=1ΣΣDTkp+kd×dtei−kd ×dt ×en−1i=1清单2：PID2的源代码。kp= 9.4514;ki = 0.69006;kd = 2.8454;invdt =5.0;dt=0.2;m= 8.0;c = 5.0;eold = 0.0;R = 0.0;s= 0.0;C1 =KP +kd *invdt;c2=kd *invdt;c3=ki *dt;while true doe=c-m;i= i+ c2 *e;R =（c1 * e）-（c3 *eold）;r= R +1;约尔德=e;M =M +0.01*r;4、如何获得更精确的PID？初始PID可以通过展开循环进行更彻底的变换。在我们的例子中，我们任意选择将其展开四次，以便每次执行时保持最后四个错误的总和。然后，我们改变操作的顺序，或者通过两两相加，或者以递增或递减的顺序。下面给出了应用于获得第三PID算法（称为PID3）的过程。让我们从展开四次积分项开始，就像静态分析器或编译器通常做的那样：in−1=in−2+ki×dt×en−1in−2=in−3+ki×dt×en−2in−3=in−4+ki×dt×en−3in−4=in−5+ki×dt×en−4我们在in中内联了前面的表达式。 我们得到：in=in−5+（ki×dt×en−4）+（ki×dt×en−3）+（ki×dt×en−2）+（ki×dt×en−1）+（ki×dt×en）（2）其中i n−5= i0+ k i× dt × n−5 e i。 等式（2）可以更简单，我们有：n−5i n= i0+ ki× dt ×ei+（ki ×dt×（e n−4 + e n−3）+e n−2）+e n−1）+e n））.i=1现在，如果我们在做了一些操作之后回到结果表达式，比如求导和因式分解，我们得到最终表达式：.1Σn−51+ k i×dt×（e n−4 + e n−3）+e n−2）+e n−1）+e n）。表示为s =。（（e n−4+ e n−3）+e n−2）+e n−1<$，k1= k p+. k d× 1dt，k2= ki× dtrn=en×+i0+ki×dt×N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）4751DT.ΣΣΣ×且k3=kd×1，则rn的最终表达式为：n−5r n= R + i + k2×s +e i，其中 R =（e n× k1）−（k3× en−1）。i=1清单3中给出了PID3的完整代码。请注意，本节中提出的转换会导致需要更多内存的更大程序虽然它允许更多的转换，但在内存很重要的某些上下文中，它可能是无关紧要的5实验结果现在让我们关注三个PID程序的执行。在我们使用GMPY2库进行多精度计算的Python实现中，我们的结果表明PID1、PID2和PID3之间存在显著差异。只要在结果的十进制值的第二位或第三位上。 为了更好地可视化这些结果，图1中给出了三种PID算法对应的曲线。我们可以观察到对应于三个PID的曲线之间的显著差异，主要在x轴的0和150之间图2显示了PID1和PID3之间的差异。这种差异是很重要的，它的计算精度很高。因此，通过使用24位、53位和50000位的尾数观察到相同的行为。PID1和PID3之间的误差在− 0之间振荡。1和0。25，而该值的范围在5和8之间。我们还观察到PID1和PID2之间的差异可以忽略不计。康-从PID1和PID3来看，我们可以注意到，代码中的一个小的语法变化确实在准确性方面产生了重要的差异。例如，让我们取PID1的以下表达式r，让我们内联三个项p，i和d，并将e因子分解。最初，r=p+i+d，我们在因式分解后得到RJ= e×。KP+ki ×dt+kd1+Idt0-我是说...K1d×dt×Eold，通过这种简单的修改，精度的差异已经很重要了，如图3所示，它给出了循环第一次迭代时r和rJ之间以及m和mJ清单3：PID3的源代码。kp= 9.4514;ki = 0.69006;kd = 2.8454;invdt =5.0;dt=0.2;R= 0.0;S = 0.0;s = 0.0;m = 8.0;c = 5.0;eold=0.0;e1 = e2 = e3 = e4 = 0.0; k1 = kp + kd * invdt; k2 = ki * dt; k3 = kd * invdt;而true确实e=c-m;i= i+ k2 *e;R =（ k1* e）-（k3* eold）;S =s+（e4+（e3+（e2+e1）;r=R+i+（k2 *S）;eold=e;e4=e3;e3=e2;e2=e1;e1=e;m=m+ 0.01 * r;52N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）47六、结论本文着重介绍了对标准PID控制器的改造。我们相信，通过使用独立于本案例研究中使用的示例程序的系统和一般规则，可以自动获得从PID1到PID2和PID3的整个转换（见第3节和第4这些规则包括表达式的内联、部分求值和循环展开。运行这三个代码时获得的结果表明，这些转换显著影响结果的准确性（几个重要因素）。目前，我们正在开发一个软件，基于前一段中提到的规则以及重写表达式的规则（结合性，交换性等）。该工具旨在将初始程序（如PID1）作为输入，并自动生成在数学上等效且更精确的其他PID程序一个关键的问题是在转换过程中发生的组合爆炸，因为在每一步都可能使用许多规则。全面的探索的等价程序集是不现实的，我们的目标是开发更易于处理的技术，使用抽象。7.576.565.55电话：+86-510 - 8888888传真：+86-510 - 8888888图1.一、三种PID算法中的测量值m测量PID 1的m1测量PID 2的m2测量PID3的m3N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）4753PID1和PID30.250.20.150.10.050-0.05电话：+86-050 -8850000传真：+86-050 - 8850000图二、PID1和PID3中测量值m之间的差异。它 RR'MM'1-71。449234-12。165627-19。748718-17。074722-16。089169-14。934349-13。892138-12。913241-12。000034-11。147227-71。449242-12。993700-21。010429-18。747597-18。077232-16。752943-15。672437-14。609431-13。609982-12。662717第七章285508第七章2855082第七章163851第七章1555713第六章966364第六章9454664第六章795617第六章7579905第六章634725第六章5772186第六章485382第六章4096887第六章346460第六章2529648第六章217328第六章1068699第六章097328五、97076910五、985856五、844142图3.第三章。 比较结果r和r′以及相应的测量值m和m′。引用[1] ANSI/IEEE。 IEEE二进制浮点运算标准，std 754-2008版，2008年。[2] 亚历山大·查普特区间斜率作为离散点变量的数值抽象域。在静态分析-第17届国际研讨会，SAS 2010，佩皮尼昂，法国，2010年9月14日至16日。第6337卷，第184-200页。Springer，2010.[3] P. Cousot和R.库索抽象解释：一个统一的格子模型，用于通过构造固定点的近似来静态分析程序。第四届ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages ， 238-252 页 ACMPress，New York，NY，1977.[4] P.lawout和R.库索通过抽象解释的程序转换框架的系统设计。《程序设计语言原理》，第178ACM，2002年。[5] D. Delmas，E. Goubault，S. Putot，J. Souyris，K. T e k kal和V. 维德琳 为了在安全关键航空电子软 件 上 实 现 FLUCTUAT 的 工 业 应 用 。 工 业 临 界 系 统 的 形 式 方 法 ， 计 算 机 科 学 讲 义 第 5825 卷 ， 第53Springer，2009年。54N. Damouche et al. /Electronic Notes in Theoretical Computer Science 317（2015）47[6] 埃里克·古博特。静态分析的抽象解释的数值程序和系统，和FLUCTUAT。在静态分析-第20届国际研讨会，SAS 2013。会议记录，第7935卷，第1-3页。Springer，2013.[7] 埃里克·古博特和西尔维·普托。有限精度计算的静态分析。验证，模型检查和抽象解释-第12届国际会议，VMCAI 2011。第6538卷，第232-247页。施普林格，2011年。[8] Arnault Ioualalen 和 Matthieu Martel 。 一 个 新 的 表 示 数 学 等 价 表 达 式 的 抽 象 域 。 在 静 态 分 析 研 讨 会（SASSpringer Verlag，2012.[9] 马蒂厄·马特尔基于语义的算术表达式转换。在Static Analysis，第14届国际研讨会，SAS 2007，Proceedings，第4634卷，第298-314页。Springer，2007年。[10] 马蒂厄·马特尔算术表达式的精确求值（特邀演讲）。电子笔记理论计算机科学，287：3[11] J. -- M. Muller，N. Brisebarre，F.deDine chin，C.- P. Jeanner od，V. Le f'evre，G. Melquiond，N. 关于我们，D. Steh l'e和S. 托雷斯手簿的浮动点A算法。BirkhauserBoston，2010.[12] Michael L.奥弗顿数值计算与IEEE浮点运算-包括一个定理，一个经验法则，和一百零一个练习。SIAM，2001年。