, 邵之江 浙江大学 控制科学与工程学院, 流程生产质量优化与控制国际联合研究中心, 杭州 310027
收稿日期：2018-09-27
基金项目：国家重点研发计划资助项目（2017YFA0700300）；NSFC-浙江两化融合联合基金资助项目（U1509209）；国家自然科学基金资助项目（61773340）；中央高校基本科研业务费专项（2018QNA5011）
作者简介：徐祖华(1976-), 男, 副教授
通信作者：赵均, 副教授, E-mail:jzhao@zju.edu.cn
摘要：针对模型预测控制（model predictive control，MPC）优化求解中占用资源较多、实时性较低且实现相对复杂的问题，该文提出了一种基于现场可编程模拟阵列（field programmable analog array，FPAA）模拟神经网络的快速模型预测控制算法。通过FPAA模拟电路来实现基于连续神经网络的二次规划求解，有效规避了离散神经网络的收敛性问题，具有求解速度快、占用资源小、简单易实现的特点；通过平移变换和尺度变换方法，解决FPAA模拟电路的信号限制。最后该文给出了FPAA模拟神经网络预测控制软硬件设计方案并通过实验验证了该算法的有效性。
关键词：快速预测控制现场可编程模拟阵列(FPAA)模拟神经网络
Fast predictive control algorithm based on an FPAA analog neural network
XU Zuhua, HUANG Yanchun, CHEN Minghao, ZHAO Jun
, SHAO ZhijiangNational Center for International Research on Quality-Targeted Process Optimization and Control, College of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China
Abstract: A fast predictive control algorithm was developed for simpler, faster MPC optimization based on an FPAA analog neural network. An FPAA analog circuit provides the quadratic programming using a continuous neural network which avoids the convergence problem of discrete neural networks in a fast, simple and flexible algorithm that uses less computational resources than previous methods. The signal constraint of the FPAA analog circuit is solved by translation and scaling. The software and hardware design of the FPAA analog neural network predictive control is presented and verified in tests that show that the algorithm is effective.
Key words: fast model predictive controlfield programmable analog array (FPAA)analog neural network
模型预测控制(model predictive control, MPC)具有处理约束和多变量耦合的优势，在工业过程中得到了广泛应用。MPC需要在线求解带约束的二次规划(quadratic programing，QP)问题，因此占用较多的计算资源和较长的计算时间，早期MPC大多采用上位机作为控制平台应用在慢采样的流程中。为了拓展MPC的应用领域，很多****致力于快速模型预测控制算法的研究。常见的几类快速模型预测控制算法研究有显式预测控制(explicit MPC, EMPC)、基于传统数值解法的MPC、基于神经网络的MPC、离散和在线相结合的MPC等。
神经网络求解二次规划与传统数值解法求解二次规划一样，都属于在线求解方法。神经网络的优势在于具有天然并行性、自适应性和电路可实现性。早期，****们[1-5]为了得到求解优化问题的神经网络模型，主要做法为引入惩罚项，将有约束优化问题等价转化为无约束优化问题，基于梯度法推导神经网络模型。因为采用了罚函数法，该类型的网络需要考虑惩罚因子的参数选择问题。这不仅使得网络实施起来更加困难，而且还只能收敛到最优解的近似值。为了得到最优解的精确值，****们[6-9]不再采用罚函数法，而是通过投影定理将不等式约束转化为等式约束，基于对偶原理和KKT(Karush-Kuhn-Tucker)条件直接构造神经网络模型。该类型神经网络可以准确求解QP问题，但是不同的网络具有不同的复杂度。目前，已经有基于神经网络的预测控制算法理论研究，但是尚不如基于传统数值解法的预测控制算法应用广泛，主要原因在于神经网络的优势需要基于模拟电路实现方式才能充分体现，而模拟电路实现神经网络困难较大。因此，目前大多将网络进行离散化，然后通过数字电路来实现。但正因为如此，神经网络相对于传统数值解法，并没有展现出明显的优势。
一般来说，模拟电路的实时性总是高于数字电路。关于模拟电路在优化领域的理论研究很早就开始了，Dennis等[10]设计了模拟电路用来求解QP问题。Chua等[11]利用模拟电路实现了Hopfield神经网络，并用于求解非线性优化问题。Vichik等[12]针对线性规划和QP的求解，设计了模拟电路。Tank和Hopfield[13]展示了模拟电路是如何用来求解优化问题的。Vichik等[14]设计了用于求解QP的模拟电路并给出了模拟电路的稳定性证明。
随着可编程模拟器件(PDA)的发展，出现了现场可编程模拟阵列(field programmable analog array，FPAA)。由于FPAA既具有模拟电路经济、简洁、低功耗以及高速的优点，又具有可编程配置能力，因此在各领域中应用逐渐增多。本文采用FPAA进行了简化对偶神经网络(simplified dual neural network, SDNN)的模拟电路设计和实现，具有较大的通用性，并进一步提出了基于FPAA模拟神经网络的快速模型预测控制算法。
1 基于状态空间的MPC算法考虑如下的线性时不变离散状态空间模型：
| $\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{x}}_{\rm{m}}}(k + 1) = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{x}}_{\rm{m}}}(k) + {\mathit{\boldsymbol{B}}_{\rm{m}}}\mathit{\boldsymbol{u}}(k),}\\{\mathit{\boldsymbol{y}}(k) = {\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{x}}_{\rm{m}}}(k).}\end{array}$ | (1) |
定义新的状态变量$\boldsymbol{x}(k)=\left[\Delta \boldsymbol{x}_{\mathrm{m}}(k)^{\mathrm{T}} \boldsymbol{y}(k)^{\mathrm{T}}\right]^{\mathrm{T}}$可以将式(1)转化为如下的状态空间增量模型为
| $\begin{array}{*{20}{c}}{\mathit{\boldsymbol{x}}(k + 1) = \mathit{\boldsymbol{Ax}}(k) + \mathit{\boldsymbol{B}}\Delta \mathit{\boldsymbol{u}}(k),}\\{\mathit{\boldsymbol{y}}(k) = \mathit{\boldsymbol{Cx}}(k).}\end{array}$ | (2) |
| $\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{Y}}_P}(k) = [\mathit{\boldsymbol{y}}{{(k + 1|k)}^{\rm{T}}}\mathit{\boldsymbol{y}}{{(k + 2|k)}^{\rm{T}}} \cdots }\\{\mathit{\boldsymbol{y}}{{(k + P|k)}^{\rm{T}}}{]^{\rm{T}}},}\end{array}$ |
| $\begin{array}{*{20}{c}}{\Delta {\mathit{\boldsymbol{U}}_M}(k) = \left[ {\Delta \mathit{\boldsymbol{u}}{{(k)}^{\rm{T}}}\Delta \mathit{\boldsymbol{u}}{{(k + 1)}^{\rm{T}}} \cdots } \right.}\\{\Delta \mathit{\boldsymbol{u}}{{(k + M - 1)}^{\rm{T}}}{]^{\rm{T}}}.}\end{array}$ |
| ${\mathit{\boldsymbol{Y}}_P}(k) = \mathit{\boldsymbol{Fx}}(k|k) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta {\mathit{\boldsymbol{U}}_M}(k).$ | (3) |
| $\boldsymbol{F}=\left[ \begin{array}{c}{C A} \\ {C A^{2}} \\ {C A^{3}} \\ {\vdots} \\ {C A^{P}}\end{array}\right],$ |
| $\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \left[ {\begin{array}{*{20}{c}}{CB}&{}&{}&{}&{}\\{CAB}&{CB}&{}&{}&{}\\{C{A^2}B}&{CAB}&{CB}&{}&{}\\ \vdots & \vdots & \vdots &{}&{}\\{C{A^{P - 1}}B}&{C{A^{P - 2}}B}&{C{A^{P - 3}}B}& \cdots &{C{A^{P - M}}B}\end{array}} \right].$ |
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{\Delta {\mathit{\boldsymbol{U}}_M}\left( k \right)} J\left( k \right) = \left\| {{\mathit{\boldsymbol{Y}}_r} - {\mathit{\boldsymbol{Y}}_P}} \right\|_{{\mathit{\boldsymbol{Q}}_y}}^2 + \left\| {\Delta {\mathit{\boldsymbol{U}}_M}\left( k \right)} \right\|_{{\mathit{\boldsymbol{Q}}_{\Delta u}}}^2,}\\{{\rm{s}}.\;{\rm{t}}.\;\;\left\{ \begin{array}{l}{\mathit{\boldsymbol{Y}}_P}\left( k \right) = \mathit{\boldsymbol{Fx}}\left( {k\left| k \right.} \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta {\mathit{\boldsymbol{U}}_M}\left( k \right),\\{\mathit{\boldsymbol{U}}_{\min }} \le {\mathit{\boldsymbol{U}}_M}\left( k \right) \le {\mathit{\boldsymbol{U}}_{\max }}.\end{array} \right.}\end{array}$ | (4) |
2 基于FPAA模拟神经网络的QP求解2.1 基于神经网络的QP求解MPC的在线优化问题式(4)可以转化为如下所示的标准QP问题：
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_{{\rm{QP}}}}} J\left( {{\mathit{\boldsymbol{x}}_{{\rm{QP}}}}} \right) = \frac{1}{2}\mathit{\boldsymbol{x}}_{{\rm{QP}}}^{\rm{T}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} + {\mathit{\boldsymbol{c}}^{\rm{T}}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}},}\\{{\rm{s}}.\;{\rm{t}}.\;\;l \le \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} \le \mathit{\boldsymbol{h}}.}\end{array}$ | (5) |
目前，MPC广泛采用的传统数值解法有内点法(interior-point method，IPM)和积极集法(active set method，ASM)。但是，传统数值解法大多涉及矩阵求逆、分解等操作，因此占用资源较多、实时性较低且实现相对复杂。神经网络具有天然并行性、自适应性和电路可实现性的优势，为实时求解大规模的QP问题提供了新的途径。本文采用由Liu等[9]提出的SDNN。
式(5)所对应的SDNN模型为
| $\left\{ {\begin{array}{*{20}{l}}{\varepsilon \frac{{{\rm{d}}\mathit{\boldsymbol{v}}}}{{{\rm{d}}t}} = - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} + \mathit{\boldsymbol{g}}\left( {\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} - \mathit{\boldsymbol{v}}} \right),}\\{{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} = \mathit{\boldsymbol{Pv}} - \mathit{\boldsymbol{s}}.}\end{array}} \right.$ | (6) |
定理2??SDNN模型的平衡输出xQP*=Pv*－ s是QP(式(5))的最优解。其中v*是SDNN模型的平衡状态点[9]。
基于定理1和定理2，可知对偶神经网络模型(式6)全局收敛至QP(式5)的最优解。因此，利用SDNN便能够求解MPC中的QP问题。
对于神经网络来求解QP问题，传统的做法是将其进行离散化，离散化方法包括Euler法、Runge-Kutta法等，然后通过数字电路来实现。但是，SDNN经过离散化后存在收敛性问题，且其收敛性受步长的影响；另外，在单核实现方式下，离散SDNN算法退变为串行计算，没有了网络的并行优势。因此，求解速度相对传统数值解法提升有限。连续神经网络没有前述问题，因此本文研究连续SDNN的模拟电路实现方式。
2.2 基于FPAA的模拟电路实现FPAA具有模拟电路经济、简洁、低功耗以及高速的优点，且支持在线编程。因此，可方便地实现相应的模拟电路功能。
FPAA基于开关电容技术，其核心是可编程模拟单元(configurable analog block，CAB)，如图 1所示。
 |
| 图 1 FPAA结构示意图 |
| 图选项 |
它的基本结构与可编程逻辑器件类似。其中，配置逻辑通过串、并行总线的在线系统编程(in-system programmability，ISP)方式接收外部输入的配置数据，然后存入存储器。输入、输出单元负责对输入、输出信号进行驱动与偏置。CAB通常由放大器和外围电容、电阻、开关等构成。其中配置数据可以控制各元件的取值及相互连接关系，从而以多种CAB的组合来实现用户需要的模拟电路功能。基于FPAA的模拟电路设计的具体过程如图 2所示。
 |
| 图 2 基于FPAA的模拟电路设计流程 |
| 图选项 |
目前，FPAA已被广泛用于模拟计算、信号调理、电路实验、仪器仪表等低成本、低功耗、小型化电子系统的设计与实现等方面。在模拟领域中展现出了巨大的潜力。
2.3 连续神经网络的FPAA实现对于MPC的优化命题式(4)，其对应的QP问题的参数为
| $\left\{ \begin{gathered} {\mathit{\boldsymbol{x}}_{{\text{QP}}}} \approx \Delta {\mathit{\boldsymbol{U}}_M}\left( k \right) \in {\mathbb{R}^{{n_\mathit{\boldsymbol{u}}}M}}, \hfill \\ \mathit{\boldsymbol{W}} = 2\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\text{T}}}{\mathit{\boldsymbol{Q}}_y}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} + {\mathit{\boldsymbol{Q}}_{\Delta u}}} \right) \in {\mathbb{R}^{{n_\mathit{\boldsymbol{u}}}M \times {n_\mathit{\boldsymbol{u}}}M}}, \hfill \\ \mathit{\boldsymbol{E}} = {\mathit{\boldsymbol{B}}_\mathit{\boldsymbol{u}}} \in {\mathbb{R}^{{n_\mathit{\boldsymbol{u}}}M \times {n_\mathit{\boldsymbol{u}}}M}}. \hfill \\ \mathit{\boldsymbol{c}} = - 2\left[ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\text{T}}}{\mathit{\boldsymbol{Q}}_y}\left( {{\mathit{\boldsymbol{Y}}_r}\left( k \right) - \mathit{\boldsymbol{Fx}}\left( {k\left| k \right.} \right)} \right)} \right] \in {\mathbb{R}^{{n_\mathit{\boldsymbol{u}}}M}}, \hfill \\ \mathit{\boldsymbol{l}} = {\mathit{\boldsymbol{U}}_{\min }} - \mathit{\boldsymbol{U}}\left( {k - 1} \right) \in {\mathbb{R}^{{n_\mathit{\boldsymbol{u}}}M}}, \hfill \\ \mathit{\boldsymbol{h}} = {\mathit{\boldsymbol{U}}_{\max }} - \mathit{\boldsymbol{U}}\left( {k - 1} \right) \in {\mathbb{R}^{{n_\mathit{\boldsymbol{u}}}M}}. \hfill \\ \end{gathered} \right.$ |
由于模拟电路本身存在工作电压范围的限制，以及不同的模拟元器件也都存在或多或少的限制。因此，需要对这些限制进行分析，结合FPAA中的模拟元器件参数的限制，可以得到如下应满足的要求：
1) 电路中传递的电信号(电压)要处于工作电压范围内：式(6)中 s、xQP、v的各元素取值范围在-3~3 V；式(6)中矩阵运算乘积(如ExQP等)的各元素取值范围在-3~3 V；|l|、|h|各元素绝对值取值范围在0.05~3 V；
2) 电路中模拟元器件参数的设置要满足一定的条件：|l|=|h|，即饱和环节g(z)关于原点呈中心对称；E、P各元素绝对值取值范围在可调节范围内。
针对上述要求，显然QP(式(5))不总是满足|l|=|h|。因此，本文首先对饱和环节g(z)进行平移变换，使得g(z)关于原点呈中心对称。
2.3.1 平移变换对于QP (式5)，令
| $\mathit{\boldsymbol{a}} = \frac{{\mathit{\boldsymbol{h}} - \mathit{\boldsymbol{l}}}}{2},\mathit{\boldsymbol{ \boldsymbol{\varDelta} }} = \frac{{\mathit{\boldsymbol{l}} + \mathit{\boldsymbol{h}}}}{2}.$ |
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_{{\rm{QP}}}}} J\left( {{\mathit{\boldsymbol{x}}_{{\rm{QP}}}}} \right) = \frac{1}{2}\mathit{\boldsymbol{x}}_{{\rm{QP}}}^{\rm{T}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} + {\mathit{\boldsymbol{c}}^{\rm{T}}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}},}\\{{\rm{s}}.\;{\rm{t}}.\;\; - \mathit{\boldsymbol{a}} \le \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} - \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} \le \mathit{\boldsymbol{a}}.}\end{array}$ | (7) |
状态方程：
| $\varepsilon \frac{{{\rm{d}}\mathit{\boldsymbol{y}}}}{{{\rm{d}}t}} = - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} + \Delta + \mathit{\boldsymbol{g}}\left( {\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{x}}_{{\rm{QP}}}} - \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} - \mathit{\boldsymbol{v}}} \right);$ | (8) |
| ${\mathit{\boldsymbol{x}}_{{\rm{QP}}}} = \mathit{\boldsymbol{Pv}} - \mathit{\boldsymbol{s}}.$ | (9) |
对于平移QP(式7)，引入尺度变换因子μ，μ是正标量。式(7)的目标函数乘以μ2，约束条件两边同乘以μ，可得如下等价的QP问题：
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_{{\rm{QP}}}}} J\left( {{\mathit{\boldsymbol{x}}_{{\rm{QP}}}}} \right) = \frac{1}{2}\mu \mathit{\boldsymbol{x}}_{{\rm{QP}}}^{\rm{T}}\mathit{\boldsymbol{W}}\mu {\mathit{\boldsymbol{x}}_{{\rm{QP}}}} + \mu {\mathit{\boldsymbol{c}}^{\rm{T}}}\mu {\mathit{\boldsymbol{x}}_{{\rm{QP}}}},}\\{{\rm{s}}.\;{\rm{t}}.\;\; - \mu \mathit{\boldsymbol{a}} \le \mathit{\boldsymbol{E}}\mu {\mathit{\boldsymbol{x}}_{{\rm{QP}}}} - \mu \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} \le \mu \mathit{\boldsymbol{a}}.}\end{array}$ | (10) |
| ${\mathit{\boldsymbol{\tilde x}}_{{\rm{QP}}}} = \mu {\mathit{\boldsymbol{x}}_{{\rm{QP}}}},\mathit{\boldsymbol{\tilde c}} = \mu \mathit{\boldsymbol{c}},\mathit{\boldsymbol{ \boldsymbol{\tilde \varDelta} }} = \mu \mathit{\boldsymbol{ \boldsymbol{\varDelta} }},\mathit{\boldsymbol{\tilde a}} = \mu \mathit{\boldsymbol{a}}.$ |
| $\begin{array}{l}\mathop {\min }\limits_{{{\mathit{\boldsymbol{\tilde x}}}_{{\rm{QP}}}}} J\left( {{{\mathit{\boldsymbol{\tilde x}}}_{{\rm{QP}}}}} \right) = \frac{1}{2}\mathit{\boldsymbol{\tilde x}}_{{\rm{QP}}}^{\rm{T}}\mathit{\boldsymbol{W}}{{\mathit{\boldsymbol{\tilde x}}}_{{\rm{QP}}}} + {{\mathit{\boldsymbol{\tilde c}}}^{\rm{T}}}{{\mathit{\boldsymbol{\tilde x}}}_{{\rm{QP}}}},\\{\rm{s}}.\;{\rm{t}}.\;\; - \mathit{\boldsymbol{a}} \le \mathit{\boldsymbol{E}}{{\mathit{\boldsymbol{\tilde x}}}_{{\rm{QP}}}} - \mathit{\boldsymbol{ \boldsymbol{\tilde \varDelta} }} \le \mathit{\boldsymbol{\tilde a}}.\end{array}$ | (11) |
状态方程：
| $\varepsilon \frac{{{\rm{d}}\mathit{\boldsymbol{\tilde v}}}}{{{\rm{d}}\mathit{\boldsymbol{t}}}} = - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{\tilde x}}_{{\rm{QP}}}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} + \mathit{\boldsymbol{g}}\left( {\mathit{\boldsymbol{E}}{{\mathit{\boldsymbol{\tilde x}}}_{{\rm{QP}}}} - \mathit{\boldsymbol{ \boldsymbol{\tilde \varDelta} }} - \mathit{\boldsymbol{\tilde v}}} \right);$ | (12) |
| ${\mathit{\boldsymbol{\tilde x}}_{{\rm{QP}}}} = \mathit{\boldsymbol{P\tilde v}} - \mathit{\boldsymbol{\tilde s}}.$ | (13) |
| ${\mathit{\boldsymbol{\tilde x}}_{{\rm{QP}}}} = \mu {\mathit{\boldsymbol{x}}_{{\rm{QP}}}},\mathit{\boldsymbol{ \boldsymbol{\tilde \varDelta} }} = \mu \mathit{\boldsymbol{ \boldsymbol{\varDelta} }},\mathit{\boldsymbol{\tilde s}} = \mu \mathit{\boldsymbol{s}},\mathit{\boldsymbol{\tilde v}} = \mu \mathit{\boldsymbol{v}}.$ |
3 基于FPAA的预测控制软硬件设计本节将FPAA模拟电路实现的连续神经网络用于MPC中，给出了FPAA模拟神经网络预测控制的软硬件设计方法。下面将分别介绍硬件部分和软件部分。
3.1 硬件部分硬件平台架构如图 3所示。采用Anadigm公司生产的FPAA以及Anadigm公司提供的AnadigmDesigner2 IDE，通过多个FPAA级联的方式，实现求解QP的模拟电路，简称模拟电路QP求解器。同时，为了实现完整的MPC仿真，数字控制部分是基于上位机实现的。由于模拟部分和数字部分之间需要进行数据传输，因此采用数据采集卡作为两者的接口电路，该数据采集卡同时具备DA/AD功能。
 |
| 图 3 硬件平台架构示意图 |
| 图选项 |
3.2 软件部分软件部分如图 4所示。MPC主程序采用Visual Studio的MFC应用程序实现。上位机和数据采集卡、上位机和FPAA开发板之间采用串口通信方式来传输数据。为了实现数据采集卡的数据传输功能，需要提供数据采集卡驱动并封装成DLL。为了实现模拟电路在线更新功能，利用AnadigmDesigner2生成C Code版本的FPAA驱动API。MATLAB提供C/C++编译器，可以将编写好的M函数编译成DLL文件供MPC主程序调用。
 |
| 图 4 软件架构示意图 |
| 图选项 |
MPC主程序采用本节设计的模拟电路来求解QP问题。基于FPAA模拟神经网络的快速模型预测控制算法详细步骤如下。
步骤1??k=0。判断是否初始化，若是，进行初始化工作，跳至步骤3；若不是，跳至步骤2。
步骤2??更新状态变量和输出变量，反馈校正。
步骤3??通过无约束优化，获得解析解ΔUM(k)：$\Delta {\mathit{\boldsymbol{U}}_M}(k) = {\left({{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_y}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}|\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} + {\mathit{\boldsymbol{Q}}_{\Delta u}}} \right)^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_y}\left({{\mathit{\boldsymbol{Y}}_r}(k) - \mathit{\boldsymbol{Fx}}(k|k)} \right)$
步骤4??判断ΔUM(k)是否满足全部的约束条件，若满足，跳至步骤8；若不满足，跳至步骤5。
步骤5??模拟电路求解QP。
a) 平移变换，得到平移SDNN；
b) 尺度变换，得到尺度SDNN；
c) 上位机通过DA模块输出模拟电压至模拟电路，以电压形式输出的参数包括 s、Δ。
步骤6??对QP决策变量$\Delta \tilde{\boldsymbol{U}}_{M}(k)$进行采样。
a) AD模块采样时刻i=1；
b) 对QP决策变量$\Delta \tilde{\boldsymbol{U}}_{M}(k)$进行AD采样，记作$\Delta \tilde{\boldsymbol{U}}_{M}(k)$，表示对$\Delta \tilde{\boldsymbol{U}}_{M}(k)$的第i次采样得到的量，采样频率为500 kHz；
c)计算相对误差：
| ${\rm{error}} = \left\{ {\begin{array}{*{20}{c}}{\frac{{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( i \right)}\left( k \right) - \Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\|}}{{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\|}},}\\{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\| \ne 0;}\\{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( i \right)}\left( k \right) - \Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\|,}\\{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\| = 0.}\end{array}} \right.$ |
| ${\rm{error}} = \left\{ {\begin{array}{*{20}{c}}{\frac{{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( i \right)}\left( k \right) - \Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\|}}{{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\|}},}\\{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\| \ne 0;}\\{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( i \right)}\left( k \right) - \Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\|,}\\{\left\| {\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( {i - 1} \right)}\left( k \right)} \right\| = 0.}\end{array}} \right.$ |
| $\Delta \mathit{\boldsymbol{\tilde U}}_M^{\left( 0 \right)}\left( k \right) = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\tilde v}}}^{\left( 0 \right)}} + \mathit{\boldsymbol{\tilde s}}.$ |
e) i=i+1，跳至b)。
步骤7??尺度反变换得到原尺度下的QP决策变量ΔUM(k)：
| $\Delta \boldsymbol{U}_{M}(k)=\frac{1}{\mu} \Delta \widetilde{\boldsymbol{U}}_{M}(k).$ |
步骤9??k=k+1，跳至步骤2。
4 实例仿真与分析为了进一步验证模拟电路QP求解器的应用效果，本文采用如下线性时不变系统[15]进行仿真。
| $\left[ {\begin{array}{*{20}{c}}{{y_1}}\\{{y_2}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{\frac{2}{{{s^2} + 2s + 2}}}&{\frac{1}{{s + 1}}}\\{\frac{3}{{{s^2} + 4s + 8}}}&{\frac{1}{{{s^2} + 2s + 5}}}\end{array}} \right]\left[ {\begin{array}{*{20}{l}}{{u_1}}\\{{u_2}}\end{array}} \right].$ |
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{\Delta {\mathit{\boldsymbol{U}}_M}(k)} J\left( {\Delta {\mathit{\boldsymbol{U}}_M}(k)} \right) = }\\{\frac{1}{2}\Delta {\mathit{\boldsymbol{U}}_M}{{(k)}^{\rm{T}}}\mathit{\boldsymbol{W}}\Delta {\mathit{\boldsymbol{U}}_M}(k) + {\mathit{\boldsymbol{c}}^{\rm{T}}}\Delta {\mathit{\boldsymbol{U}}_M}(k),}\\{{\rm{ s}}{\rm{. t}}{\rm{. }}\mathit{\boldsymbol{l}} \le \mathit{\boldsymbol{E}}\Delta {\mathit{\boldsymbol{U}}_M}(k) \le \mathit{\boldsymbol{h}}}\end{array}$ |
| $\left\{ {\begin{array}{*{20}{l}}{\varepsilon \frac{{{\rm{d}}\mathit{\boldsymbol{v}}}}{{{\rm{d}}\mathit{\boldsymbol{t}}}} = - \mathit{\boldsymbol{E}}\Delta {\mathit{\boldsymbol{U}}_M}(k) + \mathit{\boldsymbol{g}}\left( {\mathit{\boldsymbol{E}}\Delta {\mathit{\boldsymbol{U}}_M}(k) - \mathit{\boldsymbol{v}}} \right),}\\{\Delta {\mathit{\boldsymbol{U}}_M}(k) = \mathit{\boldsymbol{Pv}} - \mathit{\boldsymbol{s}}.}\end{array}} \right.$ |
| $\mathit{\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}}{0.201\;8}&{0.004\;4}&{0.001\;4}&{0.003\;8}\\{0.004\;4}&{0.226\;6}&{0.003\;4}&{0.022\;9}\\{0.001\;4}&{0.003\;4}&{0.201\;1}&{0.003\;0}\\{0.003\;8}&{0.022\;9}&{0.003\;0}&{0.219\;9}\end{array}} \right],$ |
| $\mathit{\boldsymbol{E}} = \left[ {\begin{array}{*{20}{c}}1&0&0&0\\0&1&0&0\\1&0&1&0\\0&1&0&1\end{array}} \right],$ |
| $\mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}}{4.959\;6}&{ - 0.088\;2}&{4.927\;6}&{0.165\;1}\\{ - 0.088\;2}&{4.463\;4}&{ - 0.156\;4}&{4.000\;9}\\{ - 0.032\;0}&{ - 0.068\;2}&{4.942\;5}&{ - 0.128\;0}\\{ - 0.076\;9}&{ - 0.462\;5}&{ - 0.136\;7}&{4.134\;9}\end{array}} \right].$ |
在AnadigmDesigner2中搭建电路原理图，如图 5所示。然后，利用simulator进行仿真验证，可得初始时刻对应的QP求解结果，如图 6所示。其中，4个通道(channel1、channel2、channel3、channel4)分别代表Δu1(k)、Δu1(k+1)、Δu2(k)、Δu2(k+1)。因为尺度变化因子μ=0.2，所以模拟电路QP求解器得到的最优解为(0.446 3，2.138 0，0.332 8，1.734 3)，与MATLAB中的QuadProg QP求解器得到的最优解(0.446 3，2.138 0，0.332 8，1.734 3)是一致的。
 |
| 图 5 基于FPAA的模拟电路原理图 |
| 图选项 |
 |
| 图 6 AnadigmDesigner2的simulator仿真 |
| 图选项 |
接下来分别采用本文设计的模拟电路QP求解器和QuadProg QP求解器进行MIMO系统的MPC优化控制，可得仿真结果如图 7—9所示。
 |
| 图 7 输出曲线的比较 |
| 图选项 |
 |
| 图 8 输入曲线的比较 |
| 图选项 |
 |
| 图 9 QP求解时间的比较 |
| 图选项 |
上述仿真结果显示，采用模拟电路QP求解器MPC仿真和采用QuadProg QP求解器的MPC仿真的结果几乎是一致的。平均求解时间如表 1所示。
表 1 模拟电路QP求解器与QuadProg QP求解器比较
| 求解器 | 平均求解时间 |
| 模拟电路QP求解器 | 58.36 μs |
| QuadProg QP求解器 | 3.32 ms |
表选项
采用模拟电路QP求解器与采用QuadProg QP求解器的MPC之间的误差如图 10所示。图中显示模拟电路QP求解器相对于QuadProg QP求解器，在求解QP时存在大约3.91%的误差。
 |
| 图 10 模拟电路QP求解器与QuadProg QP求解器的误差 |
| 图选项 |
模拟电路QP求解器的误差来源很多，主要可分为噪声干扰、DA/AD的模拟/数字量化误差、导线的压降、温度漂移以及电路元器件固有的误差等。这里重点考查饱和环节参数的偏差与元器件参数的舍入误差对实验结果的影响, 如图 11—12所示。
 |
| 图 11 SDNN模型参数的舍入误差对仿真结果的影响 |
| 图选项 |
 |
| 图 12 饱和环节的偏差对仿真结果的影响 |
| 图选项 |
由仿真结果可知，在模拟电路QP求解器的误差来源中，参数舍入误差、饱和环节参数偏差是影响比较大的。因此，通过提高参数的精度以及减小饱和环节的偏差，可进一步提高模拟电路QP求解器的精度。
5 结论本文提出了利用模拟神经网络求解预测控制在线过程中的QP问题。针对离散SDNN的缺点，提出了基于FPAA模拟神经网络的快速模型预测控制算法，有效规避了步长选择问题，且保证了神经网络的天然并行性。针对FPAA模拟电路信号限制，提出了平移变化方法和尺度变化方法，使得本文设计的模拟电路QP求解器的通用性显著提高。同时，给出了FPAA模拟神经网络预测控制软硬件设计方法并进行了实验验证。实验结果表明：本文提出的算法具有更快的求解速度。最后，对模拟电路进行了误差分析，指出了进一步提高模拟电路QP求解器精度的方向。
参考文献
| [1] | HOPFIELD J J. Neurons with graded response have collective computational properties like those of two-state neurons[J]. Proceedings of the National Academy of Sciences of the United States of America, 1984, 81(10): 3088-3092. DOI:10.1073/pnas.81.10.3088 |
| [2] | KENNEDY M P, CHUA L O. Neural networks for nonlinear programming[J]. IEEE Transactions on Circuits and Systems, 1988, 35(5): 554-562. DOI:10.1109/31.1783 |
| [3] | RODRIGUEZ-VAZQUEZ A, DOMINGUEZ-CASTRO R, RUEDA A, et al. Nonlinear switched capacitor "neural" networks for optimization problems[J]. IEEE Transactions on Circuits and Systems, 1990, 37(3): 384-398. DOI:10.1109/31.52732 |
| [4] | ZHANG S, CONSTANTINIDES A G. Lagrange programming neural networks[J]. IEEE Transactions on Circuits and Systems Ⅱ:Analog and Digital Signal Processing, 1992, 39(7): 441-452. DOI:10.1109/82.160169 |
| [5] | HAM F M, KOSTANIC I. Principles of neurocomputing for science and engineering[M]. New York, USA: McGraw-Hill Higher Education, 2000. |
| [6] | HU X L, ZHANG B. A new recurrent neural network for solving convex quadratic programming problems with an application to the k-winners-take-all problem[J]. IEEE Transactions on Neural Networks, 2009, 20(4): 654-664. |
| [7] | XIA Y S, FENG G, WANG J. A recurrent neural network with exponential convergence for solving convex quadratic program and related linear piecewise equations[J]. Neural Networks, 2004, 17(7): 1003-1015. DOI:10.1016/j.neunet.2004.05.006 |
| [8] | ZHANG Y N, WANG J. A dual neural network for convex quadratic programming subject to linear equality and inequality constraints[J]. Physics Letters A, 2002, 298(4): 271-278. DOI:10.1016/S0375-9601(02)00424-3 |
| [9] | LIU S B, WANG J. A simplified dual neural network for quadratic programming with its KWTA application[J]. IEEE Transactions on Neural Networks, 2006, 17(6): 1500-1510. DOI:10.1109/TNN.2006.881046 |
| [10] | DENNIS J B. Mathematical programming and electrical networks[D]. Cambridge, USA: Massachusetts Institute of Technology, 1959. |
| [11] | CHUA L, LIN G N. Nonlinear programming without computation[J]. IEEE Transactions on Circuits and Systems, 1984, 31(2): 182-188. DOI:10.1109/TCS.1984.1085482 |
| [12] | VICHIK S, BORRELLI F. Solving linear and quadratic programs with an analog circuit[J]. Computers & Chemical Engineering, 2014, 70: 160-171. |
| [13] | TANK D W, HOPFIELD J J. Simple "neural" optimization networks:An A/D converter, signal decision circuit, and a linear programming circuit[J]. IEEE Transactions on Circuits and Systems, 1986, 33(5): 533-541. DOI:10.1109/TCS.1986.1085953 |
| [14] | VICHIK S, ARCAK M, BORRELLI F. Stability of an analog optimization circuit for quadratic programming[J]. Systems & Control Letters, 2016, 88: 68-74. |
| [15] | JIAN L, LI D W, XI Y G. Implementation of dynamic matrix control on FPGA[C]//Proceedings of the 29th Chinese Control Conference (CCC). Beijing, China: IEEE, 2010: 5970-5974. |
