最优化方法大作业

发动机空燃比控制器

引言：我主要从事自动化相关研究。这里介绍我曾经接触过的发动机空燃比控制器设计中的优化问题。

发动机空燃比控制器设计中的最优化问题

AFR =a

f

m m && (1)

空燃比由方程（1）定义，在发动机运行过程中如果控制AFR 稳定在14.7可以获

得最好的动力性能和排放性能。如果假设进入气缸的空气流量a

m &可以由相关单元检测得到，则可以通过控制进入气缸的燃油流量f m &来实现空燃比的精确控制。由于实际发动机的燃油喷嘴并不是直接对气缸喷燃油，而是通过进气歧管喷燃油，这么做会在进

气歧管壁上液化形成油膜，因此不仅是喷嘴喷出的未液化部分燃油会进入气缸，油膜

蒸发部分燃油也会进入气缸，如方程（2）。这样如何更好的喷射燃油成为了一个问题。

1110101122211ττττ??

??

-??

??????????=+????????-????????????-????

?

???

??

?????????

?f f

f 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 &、fv

m &表示为液化部分燃油、fi

m &表示喷嘴喷射的燃油，在τf 、τv 、X 都已知的情况下，由现代控制理论知识，根据系统的增广状态空间模型方程（3）

0000001

1

011011114.70ττττ????-??

??????????=-+-???????????????

???????????????

??

??=??????

f

f v v a X X u +q q m y q x x x &&& (3)

其中()0

14.7?t

a

q =

y -m

&。由极点配置方法，只要设计控制器方程（4），就可以

使得y 无差的跟踪阶跃输入，那么y 也能较好的跟踪AFR *a

m /&。 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??????=---*-??????

ff

fi a f fv m m m

m m &&&&& （5） 画出最后的仿真跟踪效果图，如图1，可见坐标轮换法很好的完成了任务。

图1 f

m &对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();

P

J = [-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 ~= 0

k1 = k1 + 1;

if f1 < f2

n = 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);

else

m = x1;

x1 = x2;

f1 = f2;

x2 = m + 0.618*(n-m); f2 = wch(x2,b);

end

end

a = (m+n)/2;flag = 0;

else

k2 = k2 +1;

if F1 < F2

N = 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);

else

M = y1;

y1 = y2;

F1 = F2;

y2 = M + 0.618*(N-M);

F2 = wch(a,y2);

end

end

b = (M+N)/2;flag = 1;

end

end

k1

k2

R = [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;

end

f = k;

end