发动机空燃比控制器引言：我主要从事自动化相关研究。

这里介绍我曾经接触过的发动机空燃比控制器设计中的优化问题。

发动机空燃比控制器设计中的最优化问题AFR =afm m && (1)空燃比由方程（1）定义，在发动机运行过程中如果控制AFR 稳定在14.7可以获得最好的动力性能和排放性能。

如果假设进入气缸的空气流量am &可以由相关单元检测得到，则可以通过控制进入气缸的燃油流量f m &来实现空燃比的精确控制。

由于实际发动机的燃油喷嘴并不是直接对气缸喷燃油，而是通过进气歧管喷燃油，这么做会在进气歧管壁上液化形成油膜，因此不仅是喷嘴喷出的未液化部分燃油会进入气缸，油膜蒸发部分燃油也会进入气缸，如方程（2）。

这样如何更好的喷射燃油成为了一个问题。

1110101122211ττττ⎡⎤⎡⎤-⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤⎢⎥⎣⎦⎢⎥⎣⎦f ff v X x x u x x X x y =x && (2)其中12、,==ff fv x m x m &&=f y m &,=fi u m &这里面，表示油膜蒸发量ff m &、fvm &表示为液化部分燃油、fim &表示喷嘴喷射的燃油，在τf 、τv 、X 都已知的情况下，由现代控制理论知识，根据系统的增广状态空间模型方程（3）00000011011011114.70ττττ⎡⎤⎡⎤-⎡⎤⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥=-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤=⎢⎥⎣⎦⎣⎦ff v v a X X u +q q m y q x x x &&& (3)其中()014.7⎰taq =y -m&。

由极点配置方法，只要设计控制器方程（4），就可以使得y 无差的跟踪阶跃输入，那么y 也能较好的跟踪AFR *am /&。

12-- u =K q K x (4)这里面的12、K K 确定，可由主导极点概念降维成两个参数12C ,C ，虽然都是最终稳态无差，但是目标是使得瞬态过程中y 和阶跃输入y r 的差异尽可能的小。

所以原问题转化成了一个关于参数的12C ,C 的，以()()2=-⎰W r c y y 为目标函数的优化问题。

这里面由于无法写出参数与目标函数的显示表达式，从而也就无法使用和导数有关的优化方法。

所以我选择用坐标轮换法，由于坐标轮换法，局部寻优性能好，且不需要事先知道函数的导数。

在matlab7.10.0上对其进行仿真试验（程序见附录）。

假设进入气缸的空气质量流量已知，并且已经适当量化，假设状态值已经由状态观测器测量得到，在τf = 1s 、τv = 60ms 、X = 0.5时按上述方法设计控制器。

坐标轮换法的精度10.051ε⎡⎤=⎣⎦,a ,b其中，11001010，''⎡⎤⎡⎤==⎣⎦⎣⎦a ,b ,为12'⎡⎤⎣⎦c ,c 的初次搜索区间；基于阶跃响应的系统性能指标()()2=-⎰W r c y y ，通过坐标轮换法得到最优参数1 4.0542,0.9544,2**=-=c c 得到最优极点1 4.0542+0.9544, 4.05420.92λλ**=-=--j 544j ，将最优极点带入爱克曼公式得到反馈增益矩阵9.7820,4.4379,43.3499⎡⎤-⎣⎦=K ，进一步得到控制器：()9.7820,4.437943.349914.7⎡⎤⎡⎤⎢⎥=---*-⎣⎦⎢⎥⎣⎦fffi a f fv m m mm m &&&&& （5） 画出最后的仿真跟踪效果图，如图1，可见坐标轮换法很好的完成了任务。

图1 fm &对14.7a m &的跟踪效果从而使得发动机空燃比的变化如图2。

图2 空燃比变化time(sec)time(sec)空燃比主函数tf = 1;tv = 0.06; X=0.5;k0 = 10000000;k = 0;A = [-1/tf 0;0 -1/tv];tf0 = 1;tv0 = 0.06;X0 = 0.5;A1 = [-1/tf0 0;0 -1/tv0];B1 = [X0/tf0 (1-X0)/tv0]';C = [1 1];D = 0;AA = [A1 zeros(2,1);C 0];BB = [B1 ;0];P = zblh();PJ = [-10*sqrt(P(1)*P(1)+P(2)*P(2)) (-P(1)+P(2)*1i) (-P(1)-P(2)*1i)];K = acker(AA, BB, J)K1 = [K(1) K(2)];K2 = K(3);-10*sqrt(P(1)*P(1)+P(2)*P(2))Ac = [(A-B*K1) -B*K2;C 0];Bc = [0;0;-1];Cc = [C 0];Dc = 0;hold on;sys = ss(Ac,Bc,Cc,Dc);a1 = 0:0.01:3;b1 = 1+0*a1;a2 = 3.01:0.01:4;b2 = 1+(a2-3);a3 = 4.01:0.01:6;b3 = 2+0*a3;a4 = 6.01:0.01:10;b4 = 2+(a4-6)/2;a5 = 10.01:0.01:14;b5 = 4+0*a5;a6 = 14.01:0.01:17;b6 = 4-(a6-14);a7 = 17.01:0.01:20;b7 = 1+0*a7;t0 = [a1 a2 a3 a4 a5 a6 a7];u0 = [b1 b2 b3 b4 b5 b6 b7];x0 = [0.5 0.5 -0.1769]';[y0,t0,x] = lsim(sys,u0,t0,x0);x1 = [1 0 0 ]*x';x2 = [0 1 0]*x';x3 = [0 0 1]*x';x = [x1;x2];u = -K1*x-K2*x3;l = u0./y0';hold on;plot(t0,u0,'-');plot(t0,y0,':');xlabel('time(sec)')ylabel('Output')zblu（）坐标轮换法函数function [ R ] = zblh( )m = 0;n = 10;x2 = m + 0.618*(n-m);f2 = wch(x2,0.5);x1 = m + 0.382*(n-m);f1 = wch(x1,0.5);M = 0;N = 10;y2 = M + 0.618*(N-M);F2 = wch(1,y2);y1 = M + 0.382*(N-M);F1 = wch(1,y1);flag = 1;a = 0.5;b = 0.5;k1 = 0;k2 = 0;while(abs(abs(m-n)+abs(M-N))>0.05) if flag ~= 0k1 = k1 + 1;if f1 < f2n = x2;x2 = x1;f2 = f1;x1 = m + 0.382*(n-m);f1 = wch(x1,b);else if f1 == f2;m = x1;n = x2;x2 = m + 0.618*(n-m); f2 = wch(x2,b);x1 = m + 0.382*(n-m); f1 = wch(x1,b);elsem = x1;x1 = x2;f1 = f2;x2 = m + 0.618*(n-m); f2 = wch(x2,b);endenda = (m+n)/2;flag = 0;elsek2 = k2 +1;if F1 < F2N = y2;y2 = y1;F2 = F1;y1 = M + 0.382*(N-M);F1 = wch(a,y1);else if F1 == F2;M = y1;N = y2;y2 = M + 0.618*(N-M); F2 = wch(a,y2);y1 = M + 0.382*(N-M);F1 = wch(a,y1);elseM = y1;y1 = y2;F1 = F2;y2 = M + 0.618*(N-M);F2 = wch(a,y2);endendb = (M+N)/2;flag = 1;endendk1k2R = [a b]end目标函数function [ f ] = wch(a,b)tf = 1;tv = 0.06;X=0.5;k0 = 10000;k = 0;A = [-1/tf 0;0 -1/tv];tf0 = 1;A1 = [-1/tf0 0;0 -1/tv];B = [X/tf (1-X)/tv]';B1 = [X/tf0 (1-X)/tv]';C = [1 1];D = 0;AA = [A1 zeros(2,1);C 0];BB = [B1 ;0];J = [-10*sqrt(a*a+b*b) (-a+b*1i) (-a-b*1i)]; K = acker(AA, BB, J);K1 = [K(1) K(2)];K2 = K(3);Ac = [(A-B*K1) -B*K2;C 0];Bc = [0;0;-1];Cc = [C 0];Dc = 0;t0 = 0:0.01:4;yn = heaviside(t0);sys = ss(Ac,Bc,Cc,Dc);u0 = heaviside(t0);y0 = lsim(sys,u0,t0);for t0 = 1:400;kn = (yn(t0) - y0(t0))*(yn(t0) - y0(t0)); k = k+kn;endf = k;end。