超级会员免费看
化工优化计算中的几种方法及应用
在化工计算中，优化问题是一个常见且重要的研究领域。本文将介绍几种用于解决约束最小化问题的方法，包括相关的理论基础、MATLAB 实现以及具体的应用示例。
1. 初始问题及求解
首先考虑一个具有约束条件的最小化问题，其约束条件可重新整理如下：
- x1 ≤ -3
-10x1 + x2 ≤ -12
- x1 ≤ 40
- x2 ≤ 40
x1 ≤ 40
x2 ≤ 40
对应的矩阵
A
和向量
b
为：
A = [-1 0; -10 1; -1 0; 0 -1; 1 0; 0 1];
b = [-3 -12 40 40 40 40]';
目标函数和其梯度定义如下：
fun = @(x) 0.02*x(1)^2 + 1.2*x(2)^2 - 80;
delfun = @(x) [0.04*x(1); 2.4*x(2)];
设置相关参数：
ne = 0;
m = 6;
crit = 1e-4;
x0 = [4 5];
使用
rosencopt
函数求解：
[xopt, fopt, iter] = rosencopt(fun, delfun, x0, A, b, ne, m, crit);
求解结果为：
xopt =
3.0000 0.0000
fopt =
-79.8200
iter =
4
2. Zoutendijk’s 可行方向法
2.1 问题描述
考虑如下形式的非线性约束最小化问题：
Minimize f(x)
Subject to gi(x) ≤ 0 (i = 1, …, m)
其中
x = [x1, x2, …, xn]T
是
n
个实值变量的列向量。
2.2 关键步骤
-
定义活动集
H：
H = {j | gj(xk) + ε ≥ 0, j = 1, …, m}
-
引入人工变量
α：
α = max{∇f(xk)T d, ∇gj(xk)T d, j ∈ H}
- 求解子问题 ：
Minimize α
Subject to ∇f(xk)T d ≤ α, ∇gj(xk)T d ≤ α, j ∈ H, -1 ≤ di ≤ 1 (i = 1, …, n)
- 步长问题 ：
Minimize f(α) = f(xk + αd)
Subject to gj(α) = gj(xk + αd) ≤ 0 (i = 1, …, m)
可重写为：
Minimize -β = α
Subject to ∇f(xk)T s + β ≤ c0, ∇gj(xk)T s + θjβ ≤ cj (j ∈ H), 0 ≤ si ≤ 2 (i = 1, …, n), β ≥ 0
其中
c0
和
cj
为相关系数。
2.3 算法步骤总结
-
选择一个满足约束条件的可行起始点
x0。 -
确定活动集
H。 -
求解最小化问题（最小化
-β = α）以找到β*和d*。 -
如果
β* = 0，停止算法。 -
否则，找到最优步长
αk并更新当前点为xk+1 = xk + αkd*，然后返回步骤 2。
2.4 MATLAB 实现
使用
zoutopt
函数实现 Zoutendijk’s 可行方向法，其基本语法为：
[xopt, fopt, iter] = zoutopt(funz, delf, delg, x0, xl, xu, nc, ncs, crit, kmax);
其中各参数含义如下：
| 参数 | 含义 |
| ---- | ---- |
|
funz
| 目标函数和约束条件 |
|
delf
| 目标函数的梯度 |
|
delg
| 约束条件的梯度 |
|
x0
| 起始点向量 |
|
xl
|
x
的下限 |
|
xu
|
x
的上限 |
|
nc
| 活动约束的数量 |
|
ncs
| 约束的总数 |
|
crit
| 停止准则 |
|
kmax
| 最大可能迭代次数 |
|
xopt
| 最优解 |
|
fopt
| 最优解处的函数值 |
|
iter
| 迭代次数 |
以下是
zoutopt
函数的具体实现：
function [xopt, fopt, iter] = zoutopt(funz, delf, delg, x0, xl, xu, nc, ncs, crit, kmax)
mcrit = 2e-3;
gcrit = 1e-4;
maxac = 20;
maxs = 30;
nv = length(x0);
x = x0;
fg = funz(x);
f = fg(1);
g = fg(2);
ic = 0;
fold = f;
iter = 0;
f0 = f;
while (1)
iter = iter + 1;
[nc, na] = actcont(ncs, crit, g); % 活动约束
if (iter > kmax)
disp('Maximim possible iterations exceeded.');
break;
end
[d, dn, beta] = dirvec(delf, delg, crit, nv, x, nc, xl, xu, mcrit); % 方向向量
if (abs(dn) < crit | abs(beta) < crit)
break;
end
d = d/dn;
alpha = linsr(funz, delf, nv, ncs, x, nc, na, xl, xu, d, x0, maxs, gcrit, crit); % 计算步长
x = x0 + alpha*d;
fg = funz(x);
f = fg(1);
x0 = x;
if (abs(f - fold) < crit)
ic = ic + 1;
if (ic == 2)
break;
end
else
ic = 0;
end
fold = f;
end
fg = funz(x);
f = fg(1);
g = fg(2);
[nc, na] = actcont(ncs, crit, g);
xopt = x;
fopt = f;
end
function [nc, na] = actcont(ncs, mcrit, g)
for i = 1:ncs
na(i) = i;
end
nc = 0;
for k = 1:ncs
if (g(k) > -mcrit)
nc = nc + 1;
ntemp = na(k);
na(k) = na(nc);
na(nc) = ntemp;
end
end
end
function [d, dn, beta] = dirvec(delf, delg, crit, nv, x, nc, xl, xu, mcrit)
df = delf(x)'; % 行向量
if (nc > 0)
A = delg(x)';
end
df = df/norm(df); % 归一化
if (nc > 0)
for j = 1:nc
A(j,:) = A(j,:)/sqrt(A(j,:)*A(j,:)');
end
end
% 活动边界
for k = 1:nv
if (xl(k) - x(k) + mcrit >= 0)
nc = nc + 1;
A(nc,1:nv) = 0;
A(nc,k) = -1;
end
end
for k = 1:nv
if (x(k) - xu(k) + mcrit >= 0)
nc = nc + 1;
A(nc,1:nv) = 0;
A(nc,k) = 1;
end
end
if (nc == 0)
beta = 1;
d = -df;
dn = norm(d);
return;
end
[beta, d] = simpx(nc, nv, df, A);
dn = sqrt(d*d'); % 可行方向
end
function alpha = linsr(funz, delf, nv, ncs, x, nc, na, xl, xu, d, x0, maxs, gcrit, crit)
nlarge = 1e40;
c = max(abs(xu - xl));
for k = 1:nv
if (abs(d(k))*nlarge > c)
if (d(k) < 0)
cn = (xl(k) - x(k)) / d(k);
if (cn < nlarge)
nlarge = cn;
end
else
cn = (xu(k) - x(k)) / d(k);
if (cn < nlarge)
nlarge = cn;
end
end
end
end
abet = nlarge;
x = x0 + abet * d;
fg = funz(x);
f = fg(1);
g = fg(2);
gmax = max(g);
inda = 1;
if (gmax <= 0)
inda = 0;
end
if (inda == 0)
amax = abet;
else
x1 = 0;
[xm, fm] = nears(funz, x1, abet, x, d, x0, maxs, crit);
amax = xm;
end
x = x0 + amax*d;
df = delf(x)';
sdr = df*d';
if (sdr <= 0)
alpha = amax;
return;
end
a1 = 0;
a2 = amax;
adif = a2 - a1;
while ((a2 - a1) > crit * adif)
am = (a1 + a2)/2;
x = x0 + am*d;
df = delf(x)';
sdr = df*d';
if (sdr == 0)
break;
end
if (sdr < 0)
a1 = am;
elseif (sdr > 0)
a2 = am;
end
end
alpha = a1;
end
function [xm, fm] = nears(funz, xa, xb, x, d, x0, maxs, crit)
miter = 0;
while (1)
xm = (xa + xb)/2;
miter = miter + 1;
if (miter > maxs)
xm = xa;
fm = fa;
return;
end
x = x0 + xm*d;
fg = funz(x);
f = fg(1);
g = fg(2);
gmax = max(g);
fm = f;
if (gmax <= 0 & gmax >= -crit)
return;
end
if (gmax < 0)
xa = xm;
fa = fm;
else
xb = xm;
end
end
end
function [beta, d] = simpx(nc, nv, df, A)
Bg = 1e2;
nrow = nc + nv + 2;
nm = nc + nv + 1;
Bm(1:nrow) = 0;
Bm(1) = sum(df);
for j = 1:nc
Bm(j+1) = 0;
for k = 1: nv
Bm(j+1) = Bm(j+1) + A(j,k);
end
end
for k = nc + 2:nrow - 1
Bm(k) = 2;
end
for k = 1: nm
Bs(k) = nv + k + 1;
end
ncol = nv + nm + 1;
for k = 1:nm
if (Bm(k) < 0)
ncol = ncol + 1;
Bs(k) = -ncol;
end
end
Am(1:nrow,1:ncol) = 0;
Am(1,1:nv) = df;
Am(1,nv+1) = 1;
for k = 1:nc
for j = 1: nv
Am(k+1,j) = A(k,j);
end
Am(k+1,nv+1) = 1;
end
mi = 0;
for k = nc+2: nrow-1
mi = mi + 1;
Am(k,mi) = 1;
end
Am(nrow,nv+1) = -1;
for k = 1: nm
Am(k,nv+k+1) = 1;
end
nt = nv + nm + 1;
for k = 1: nm
if (Bm(k) < 0)
nt = nt + 1;
Bm(k) = -Bm(k);
for j = 1:ncol
Am(k,j) = -Am(k,j);
end
Am(k,nt) = 1;
Am(nrow,nt) = Bg;
end
end
for k = 1:nm
if (Bs(k) < 0)
Bs(k) = -Bs(k);
for j = 1: ncol
Am(nrow,j) = Am(nrow,j) - Bg*Am(k,j);
end
Bm(nrow) = Bm(nrow) - Bg*Bm(k);
end
end
while (1)
nf0 = 0;
for k = 1:ncol
if (Am(nrow, k) < 0)
nf0 = 1;
break;
end
end
if (nf0 == 0)
break;
end
c = Bg;
for j = 1:ncol
if (Am(nrow, j) < c)
c = Am(nrow, j);
iv = j;
end
end
ik = 0;
jk = 0;
for k = 1:nrow - 1
if (Am(k,iv) > 0)
jk = jk + 1;
c1 = Bm(k)/(Am(k,iv) + 1e-10);
if (jk == 1)
c = c1;
jp = k;
else
if (c1 < c)
c = c1;
jp = k;
end
end
ik = 1;
end
end
Bs(jp) = iv;
if (ik == 0)
disp('Unbounded objective function.');
break;
end
c1 = 1/Am(jp,iv);
Bm(jp) = c1*Bm(jp);
for j = 1:ncol
Am(jp,j) = c1*Am(jp,j);
end
for k = 1: nrow
if (k ~= jp)
c2 = Am(k,iv);
for j = 1: ncol
Am(k,j) = Am(k,j) - c2*Am(jp,j);
end
Bm(k) = Bm(k) - c2*Bm(jp);
end
end
end
d(1:nv) = -1;
for k = 1: nm
for j = 1: nv
if (j == Bs(k))
d(j) = Bm(k) - 1;
end
end
end
beta = Bm(nrow);
end
2.5 示例
求解以下约束最小化问题：
Minimize f(x) = -(1.2x1 + 3x2)
Subject to g(x) = x1^2 + 6x2^2 - 1 ≤ 0, x1 ≥ 0, x2 ≥ 0, x1 ≤ 10, x2 ≤ 10
起始点为
x0 = [1, 0]
，无活动约束（
nc = 0
）。
定义目标函数和约束条件的函数：
function fg = funz(x)
fg(1) = -(1.2*x(1) + 3*x(2));
fg(2) = x(1)^2 + 6*x(2)^2 - 1;
end
定义目标函数和约束条件的梯度：
delf = @(x) [-1.2; -3];
delg = @(x) [2*x(1); 12*x(2)];
设置参数并求解：
x0 = [1 0];
xl = [0 0];
xu = [10 10];
nc = 0;
ncs = 1;
crit = 1e-4;
kmax = 1e3;
[xopt, fopt, iter] = zoutopt(@funz, delf, delg, x0, xl, xu, nc, ncs, crit, kmax);
求解结果为：
xopt =
0.6998 0.2916
fopt =
-1.7146
iter =
5
3. 广义约化梯度（GRG）法
3.1 问题描述
考虑如下形式的非线性约束最小化问题：
Minimize f(x)
Subject to hi(x) ≤ 0 (i = 1, …, m), xl ≤ x ≤ xu
通过引入非负松弛变量，可将不等式约束转化为等式约束，问题可重写为：
Minimize f(x)
Subject to gi(x) = 0 (i = 1, …, m), xl ≤ x ≤ xu
3.2 方法原理
GRG 方法基于使用等式约束消除变量。通常将
n + m
个设计变量分为两组：独立变量集
[Y]
和依赖变量集
[Z]
。
目标函数和约束函数的一阶导数为：
df(x) = ∑(i=1 to n-l) (∂f/∂yi) dyi + ∑(i=1 to m+l) (∂f/∂zi) dzi = ∇fY dY + ∇fZ dZ
dg(x) = ∑(i=1 to n-l) (∂g/∂yi) dyi + ∑(i=1 to m+l) (∂g/∂zi) dzi
当
Y
固定时，
g(x) + dg(x) = 0
，可得：
dZ = -B^(-1) {g(x) + C dY}
若约束在向量
x
处满足，即
g(x) = 0
，则
dZ = -B^(-1) C dY
。
方向向量
d = [dz, dy]T
，代入目标函数的导数可得：
df(x) = ∇fY dY + ∇fZ dZ = ∇fY dY - ∇fZ B^(-1) C dY
系数矩阵称为广义约化梯度
GR
：
GR = ∇fY - ∇fZ B^(-1) C
新点的更新公式为：
xk+1 = xk + αk d
3.3 算法步骤总结
-
选择一个满足
g(x0) = 0且xL ≤ x0 ≤ xU的起始点x0。 -
确定基本变量和非基本变量，并计算约化梯度
GR。 -
找到方向向量
d。若∥d∥ ≤ ε（ε为非常小的数），停止算法。 -
计算当前步长
αk并确定xk+1。 -
检查新点
xk+1的可行性：若该点可行，令xk = xk+1并返回步骤 2；否则，计算修正点xc,k+1。若f(xc,k+1) < f(xk)，令xk = xc,k+1并返回步骤 2；若f(xc,k+1) > f(xk)，令αk = αk/2，计算xk+1并重复步骤 5。
3.4 MATLAB 实现
使用
grgopt
函数实现广义约化梯度法，其基本语法为：
[xopt, fopt, iter] = grgopt(grgfun, delgrgf, delgrgg, x0, xl, xu, kmax, crit);
其中各参数含义如下：
| 参数 | 含义 |
| ---- | ---- |
|
grgfun
| 目标函数和约束条件 |
|
delgrgf
| 目标函数的梯度 |
|
delgrgg
| 约束条件的梯度 |
|
x0
| 起始点向量 |
|
xl
|
x
的下限 |
|
xu
|
x
的上限 |
|
kmax
| 最大可能迭代次数 |
|
crit
| 停止准则 |
|
xopt
| 最优解 |
|
fopt
| 最优解处的函数值 |
|
iter
| 迭代次数 |
以下是
grgopt
函数的具体实现：
function [xopt, fopt, iter] = grgopt(grgfun, delgrgf, delgrgg, x0, xl, xu, kmax, crit)
dcrit = crit/10;
fcrit = crit*1e-2;
x = x0;
[f g] = grgfun(x);
[ncs nv] = size(delgrgg(x));
for k = 1: ncs
if (abs(g(k)) > crit)
disp('Infeasible starting point.');
break;
end
end
indc = 0;
fold = f;
for iter = 1:kmax
df = delgrgf(x);
dg = delgrgg(x);
for k = 1:nv
nbr(k) = k;
end
[nbr,df,dg] = redgr(nbr,nv,ncs,x,xl,xu,df,dg, crit); % 基本变量
dm = 0.;
for k = ncs+1:nv
nk = nbr(k);
d(nk) = -df(nk);
if (x(nk) < xl(nk)+crit & df(nk) > 0)
d(nk) = 0;
end
if (x(nk) > xu(nk)-crit & df(nk) < 0)
d(nk) = 0;
end
if (dm < abs(d(nk)))
dm = abs(d(nk));
end
end
if (dm < crit)
break;
end
for k = 1: ncs
nk = nbr(k);
d(nk) = 0;
for j = ncs+1:nv
nj = nbr(j);
d(nk) = d(nk) - dg(k,nj)*d(nj);
end
end
x0 = x;
alpha = findstep(nv,x,xl,xu,d);
x = x0;
for k = 1: nv
xtemp(k) = x(k) + alpha*d(k);
end
[ftemp g] = grgfun(xtemp);
df = delgrgf(xtemp);
dg = delgrgg(xtemp);
fup = 0;
for k = 1:nv
fup = fup + df(k)*d(k);
end
if (fup > 0 | ftemp > f)
c = 0;
x0 = x;
[alpha,ftemp] = goldsec(nv,ncs,c,alpha,dcrit,d,x);
x = x0;
for k = 1:nv
xtemp(k) = x(k) + alpha*d(k);
end
end
[c g] = grgfun(xtemp);
ag = abs(g(1));
for k = 1:ncs
if (ag < abs(g(k)))
ag = abs(g(k));
end
end
% 若不可行，应用牛顿法
if (ag > crit)
blim = 20;
for bi0 = 1:blim
bnew = 0;
bi2 = 12;
for bi1 = 1:bi2
bnew = bnew + 1;
df = delgrgf(xtemp);
dg = delgrgg(xtemp);
[g] = rednewton (dg, g, nbr, ncs);
bf = 0;
for k = 1: ncs
nk = nbr(k);
xtemp(nk) = xtemp(nk) - g(k);
if ((xtemp(nk) < xl(nk)) | (xtemp(nk) > xu(nk)))
bf = 1;
break;
end
end
if (bf == 1)
break;
end
[c g] = grgfun(xtemp);
ag1 = abs(g(1));
for k = 1: ncs
if (ag1 < abs(g(k)))
ag1 = abs(g(k));
end
end
if (bf == 0 & ag1 < crit)
break;
end
bf = 1;
if (bnew > 3 | ag1 > ag)
break;
end
ag = ag1;
end
if (bf == 0)
[ftemp, g] = grgfun(xtemp);
if (ftemp < f)
break;
end
end
alpha = alpha/2;
for k = 1:nv
xtemp(k) = x(k) + alpha*d(k);
end
[c g] = grgfun(xtemp);
ag = abs(g(1));
for k = 1:ncs
if (ag < abs(g(k)))
ag = abs(g(k));
end
end
end
end
if (bf==0 & ftemp<f)
for k = 1: nv
x(k) = xtemp(k);
end
f = ftemp;
end
fcrit = abs(f)*fcrit + fcrit;
if (abs(f - fold) < fcrit)
indc = indc + 1;
if (indc > 1)
break;
end
else
indc = 0;
end
fold = f;
end
xopt = x;
fopt = f;
end
%------------------------------------------------------------
function [x4,ft] = goldsec(nv,ncs,x1,x4,dcrit,d,x)
tau = (sqrt(5)-1)/2;
x2 = tau*x1 + (1-tau)*x4;
for k = 1:nv
xtemp(k) = x(k) + x2*d(k);
end
[f2, g] = grgfun(xtemp);
for count = 1:100
x3 = tau*x4 + (1-tau)*x1;
for k = 1:nv
xtemp(k) = x(k) + x3*d(k);
end
[f3, g] = grgfun(xtemp);
if (f2 < f3)
x4 = x1;
x1 = x3;
else
x1 = x2;
x2 = x3;
f2 = f3;
end
if (abs(x4 - x1) <= dcrit)
break;
end
end
x4 = x2;
ft = f2;
end
%------------------------------------------------------------
function [g] = rednewton(dgr, g, nbr, ncs)
for k = 1: ncs-1
nk = nbr(k);
for jk = k+1:ncs
c = dgr(jk,nk)/dgr(k,nk);
for j = k + 1: ncs
nj = nbr(j);
dgr(jk,nj) = dgr(jk,nj) - c*dgr(k,nj);
end
g(jk) = g(jk) - c*g(k);
end
end
g(ncs) = g(ncs)/dgr(ncs,nbr(ncs));
for jm = 1:ncs-1
jk = ncs - jm;
im = nbr(jk);
c = 1/dgr(jk,im);
g(jk) = c*g(jk);
for k = jk + 1: ncs
nk = nbr(k);
g(jk) = g(jk) - c*dgr(jk,nk)*g(k);
end
end
end
%------------------------------------------------------------
function [nbr,df,dg] = redgr(nbr,nv,ncs,x,xl,xu,df,dg,xcrit)
for k = 1: ncs
nk = nbr(k);
kcount = 0;
for j = k: nv
nj = nbr(j);
if (x(nj) > xl(nj)+xcrit & x(nj) < xu(nj)-xcrit)
kcount = kcount + 1;
if (kcount == 1)
pivot = dg(k,nj);
jpv = j;
else
if (abs(pivot) < abs(dg(k,nj)))
pivot = dg(k,nj);
jpv = j;
end
end
end
end
nbr(k) = nbr(jpv);
nbr(jpv) = nk;
nk = nbr(k);
if (k == ncs)
break;
end
for jk = k+1: ncs
dgratio = dg(jk,nk)/dg(k,nk);
for j = k+1:nv
nj = nbr(j);
dg(jk,nj) = dg(jk,nj) - dgratio*dg(k,nj);
end
end
end
nbs = nbr(ncs);
for j = ncs + 1: nv
nj = nbr(j);
dg(ncs,nj) = dg(ncs,nj)/dg(ncs,nbs);
for jm = 1:ncs-1
jk = ncs-jm;
km = nbr(jk);
dgratio = 1/dg(jk,km);
dg(jk,nj) = dgratio*dg(jk, nj);
for k = jk+1: ncs
nk = nbr(k);
dg(jk,nj) = dg(jk,nj) - dgratio*dg(jk,nk)*dg(k,nj);
end
end
end
% 计算约化梯度
for jk = ncs+1:nv
km = nbr(jk);
for j = 1: ncs
nj = nbr(j);
df(km) = df(km) - dg(j,km)*df(nj);
end
end
end
%------------------------------------------------------------
function [alpha] = findstep(nv, x, xl, xu, d)
% 计算最大步长
kcount = 0;
for j = 1:nv
au = xu(j) - x(j); % 检查上界
if (d(j) > 1E-30*(au+1))
kcount = kcount + 1;
ad = au/d(j);
if (kcount == 1)
alpha = ad;
end
if (kcount > 1)
if (alpha > ad)
alpha = ad;
end
end
end
al = x(j) - xl(j); % 检查下界
if (-d(j) > 1E-30*(al + 1))
kcount = kcount + 1;
ad = -al/d(j);
if (kcount == 1)
alpha = ad;
end
if (kcount > 1)
if (alpha > ad)
alpha = ad;
end
end
end
end
end
3.5 示例
求解以下约束最小化问题：
Minimize f(x) = x1^2 + x2^2 + 2.3x3^2 - 1.2x4^2 - 4x1 - 6x2 - 20x3 + 6x4 + 100
Subject to
h1(x) = x1^2 + x2^2 + x3^2 + x4^2 + x1 - x2 + x3 - x4 ≤ 7
h2(x) = x1^2 + 2x2^2 + x3^2 + 2x4^2 - x1 - 4x2 + 2 ≤ 11
h3(x) = x1^2 + x2^2 + x3^2 + 2x1 - x2 - x4 ≤ 6
引入松弛变量将不等式约束转化为等式约束：
g1(x) = x1^2 + x2^2 + x3^2 + x4^2 + x1 - x2 + x3 - x4 + x5 - 7 = 0
g2(x) = x1^2 + 2x2^2 + x3^2 + 2x4^2 - x1 - 4x2 + 2 + x6 - 11 = 0
g3(x) = x1^2 + x2^2 + x3^2 + 2x1 - x2 - x4 + x7 - 6 = 0
起始点为
x0 = [0 0 0 0 7 11 6]
，设置
kmax = 1000
和
crit = 1×10−4
。
定义目标函数和约束条件的函数：
function [f g] = grgfun(x)
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);
x5 = x(5);
x6 = x(6);
x7 = x(7); % 松弛变量
f = x1^2 + x2^2 + 2.3*x3^2 - 1.2*x4^2 - 4*x1 - 6*x2 - 20*x3 + 6*x4 + 100;
g(1) = x1^2 + x2^2 + x3^2 + x4^2 + x1 - x2 + x3 - x4 + x5 - 7;
g(2) = x1^2 + 2*x2^2 + x3^2 + 2*x4^2 - x1 - x4 + x6 - 11;
g(3) = 2*x1^2 + x2^2 + x3^2 + 2*x1 - x2 - x4 + x7 - 6;
end
定义目标函数和约束条件的梯度：
function df = delgrgf(x)
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);
df = zeros(1,length(x));
df(1) = 2.3*x1 - 4;
df(2) = 2*x2 - 6;
df(3) = 4.6*x3 - 20;
df(4) = -2.4*x4 + 6;
end
function dg = delgrgg(x)
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);
dg = zeros(3,length(x));
dg(1,1) = 2*x1 + 1;
dg(1,2) = 2*x2 - 1;
dg(1,3) = 2*x3 + 1;
dg(1,4) = 2*x4 - 1;
dg(1,5) = 1;
dg(2,1) = 2*x1 - 1;
dg(2,2) = 4*x2 - 4;
dg(2,3) = 2*x3;
dg(2,4) = 4*x4;
dg(2,6) = 1;
dg(3,1) = 4*x1 + 2;
dg(3,2) = 2*x2 - 1;
dg(3,3) = 2*x3;
dg(3,4) = -1;
dg(3,7) = 1;
end
设置参数并求解：
x0 = [0 0 0 0 7 11 6];
xl = [0 0 0 0 0 0 0];
xu = [10 10 10 10 10 10 10];
crit = 1e-4;
kmax = 1e3;
[xopt, fopt, iter] = grgopt(@grgfun, @delgrgf, @delgrgg, x0, xl, xu, kmax, crit);
总结
本文介绍了化工计算中解决约束最小化问题的几种方法，包括初始问题的求解、Zoutendijk’s 可行方向法和广义约化梯度（GRG）法。每种方法都有其独特的原理和应用场景。
-
初始问题求解
：通过整理约束条件，使用
rosencopt函数可以直接求解简单的约束最小化问题。 - Zoutendijk’s 可行方向法 ：适用于非线性约束最小化问题，通过定义活动集、引入人工变量和求解子问题来找到可行方向和步长，最终得到最优解。
- 广义约化梯度（GRG）法 ：主要处理非线性等式约束，通过消除变量和计算约化梯度来更新点，逐步逼近最优解。
这些方法在 MATLAB 中都有相应的实现，通过具体的示例展示了如何使用这些方法解决实际问题。在实际应用中，可以根据问题的特点选择合适的方法进行求解。
流程图
以下是 Zoutendijk’s 可行方向法和广义约化梯度（GRG）法的流程图：
graph TD;
A[选择可行起始点x0] --> B[确定活动集H];
B --> C[求解子问题得β*和d*];
C --> D{β* = 0?};
D -- 是 --> E[停止算法];
D -- 否 --> F[找到最优步长αk];
F --> G[更新当前点xk+1 = xk + αkd*];
G --> B;
graph TD;
A[选择起始点x0] --> B[确定基本和非基本变量，计算GR];
B --> C[找到方向向量d];
C --> D{∥d∥ ≤ ε?};
D -- 是 --> E[停止算法];
D -- 否 --> F[计算步长αk并确定xk+1];
F --> G{新点xk+1可行?};
G -- 是 --> H[xk = xk+1];
H --> B;
G -- 否 --> I[计算修正点xc,k+1];
I --> J{f(xc,k+1) < f(xk)?};
J -- 是 --> K[xk = xc,k+1];
K --> B;
J -- 否 --> L[αk = αk/2，计算xk+1];
L --> G;
这些流程图清晰地展示了两种算法的步骤和决策过程，有助于理解算法的执行逻辑。
扫一扫
766
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
