一.線性規(guī)劃
1.算法
【1】通用代碼
[x,y]=linprog(c,A,b,Aeq,beq,lb,ub)
例如：
max需要加負(fù)號(hào)變成min、>=需要加負(fù)號(hào)變成<=
matlab
（1）基于求解器
%? ? LP1_1_1.m
f=[-2;-3;5];
a=[-2,5,-1;1,3,1];b=[-10;12];
aeq=[1,1,1];beq=7;
[x,y]=linprog(f,a,b,aeq,beq,zeros(3,1));
x,y=-y? ? %最大值
x =
6.4286
0.5714
0
y =
14.5714
% 說明：在x1 x2 x3 = 6.4286 0.5714 0 的情況下，maxz = 14.5714（2）基于問題
con中根據(jù)符號(hào)分類
% LP1_1_2.m
clc,clear
prob=optimproblem('ObjectiveSense','max')
x=optimvar('x',3,'LowerBound',0);
prob.Objective=2*x(1)+3*x(2)-5*x(3);
prob.Constraints.con1=x(1)+x(2)+x(3)==7;
prob.Constraints.con2=2*x(1)-5*x(2)+x(3)>=10;
prob.Constraints.con3=x(1)+3*x(2)+x(3)<=12;
[sol,fval,flag,out]=solve(prob),sol,x;
x=sol.x,y=fvalpython
#? ? 線性規(guī)劃模型.py
from scipy.optimize import linprog
import numpy as np
c=np.array([-2,-3,5])
a=[[-2,5,-1],[1,3,1]]
b=[[-10],[12]]
aeq=[[1,1,1]]
beq=[7]
lb=[0,0,0]
ub=[None]*len(c)
bound=tuple(zip(lb,ub))
res=linprog(c,a,b,aeq,beq,bound,method='simplex',options={"disp":True})
print("最優(yōu)解：",res.x)
print("目標(biāo)函數(shù)最小值：",res.fun)
print("目標(biāo)函數(shù)最大值：",-res.fun)
print(res)
最優(yōu)解： [6.42857143 0.57142857 0. ]
目標(biāo)函數(shù)最小值： -14.571428571428571
目標(biāo)函數(shù)最大值： 14.571428571428571
con: array([0.])
fun: -14.571428571428571
message: 'Optimization terminated successfully.'
nit: 3
slack: array([0. , 3.85714286])
status: 0
success: True
x: array([6.42857143, 0.57142857, 0. ])【2】可以轉(zhuǎn)化為線性規(guī)劃
（1）絕對(duì)值
%? ? LP_abs.m
clc,clear
c=[1,2,3,4]';
a=[1,-1,-1,1;1,-1,1,-3;1,-1,-2,3];
b=[-2,-1,-1/2]';
prob=optimproblem;
u=optimvar('u',4,'LowerBound',0);
v=optimvar('v',4,'LowerBound',0);
prob.Objective=sum(c'*(u+v));
prob.Constraints.con=a*(u-v)<=b;
[sol,fval,flag,out]=solve(prob)
x=sol.u-sol.v,minz=fval
x =
-2
0
0
0
minz =
2.0000（2）min(max(q*x))
（見風(fēng)投案例模型二）
2.風(fēng)投案例
【0】題目描述
【1】模型一
模型一：設(shè)定風(fēng)險(xiǎn)度的最大接受值，在不太冒險(xiǎn)的情況下選擇最大收益。
（1）matlab
%? ? LP_fengtou_1.m
clc,clear
t=0;hold on
while t<0.05
c=[-0.05,-0.27,-0.19,-0.185,-0.185];
A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])];
b=t*ones(4,1);
Aeq=[1,1.01,1.02,1.045,1.065];beq=1;
lb=zeros(5,1);
[x,Q]=linprog(c,A,b,Aeq,beq,lb);
Q=-Q;
plot(t,Q,'*k');
t=t+0.001;
end
xlabel('t'),ylabel('Q')
選擇曲線的轉(zhuǎn)折點(diǎn)作為最優(yōu)投資組合。
輸出轉(zhuǎn)折點(diǎn)見模型一python代碼：
鼠標(biāo)在圖上確定大致坐標(biāo)，在循環(huán)中加入條件輸出。
（2）python
#? ? 風(fēng)投案例.py
from scipy.optimize import linprog
import numpy as np
from numpy import ones,diag,c_,zeros
import matplotlib.pyplot as plt
plt.rc('font',size=16)
c=np.array([-0.05,-0.27,-0.19,-0.185,-0.185])
A=c_[zeros(4),diag([0.025,0.015,0.055,0.026])]
Aeq=[[1,1.01,1.02,1.045,1.065]]
beq=[1]
a=0
aa=[]
ss=[]
while a<0.05:
b=ones(4)*a
res = linprog(c, A, b, Aeq, beq)
x=res.x
Q=-res.fun
aa.append(a)
ss.append(Q)
#輸出轉(zhuǎn)折點(diǎn)
if a==0.006:
print("x=", x)
print("Q=", Q)
a=a+0.001
plt.plot(aa,ss,'r*')
plt.xlabel('$a$')
plt.ylabel('$Q$',rotation=90)
plt.savefig('模型一.png',dpi=500)
plt.show()
x= [0. 0.24 0.4 0.10909091 0.22122066]
Q= 0.20190763977806234（萬元）
【2】模型二
模型二的目標(biāo)函數(shù)：使各項(xiàng)目中風(fēng)險(xiǎn)最大的一項(xiàng)風(fēng)險(xiǎn)最小。
import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp
plt.rc('font',size=16)
x=cp.Variable(6,pos=True)
obj=cp.Minimize(x[5])
a1=np.array([0.025,0.015,0.055,0.026])
a2=np.array([0.05,0.27,0.19,0.185,0.185])
a3=np.array([1,1.01,1.02,1.045,1.065])
k=0.05
kk=[]
ss=[]
while k<0.27:
con=[cp.multiply(a1,x[1:5])-x[5]<=0,
0.05*x[0]+0.27*x[1]+0.19*x[2]+0.185*x[3]+0.185*x[4]>=k,
x[0]+1.01*x[1]+1.02*x[2]+1.045*x[3]+1.065*x[4]==1]
prob=cp.Problem(obj,con)
prob.solve(solver='CVXOPT')
#輸出轉(zhuǎn)折點(diǎn)
if abs(k-0.21)<0.0005:
print("x=",x.value)
print("風(fēng)險(xiǎn)=",prob.value)
kk.append(k)
ss.append(prob.value)
k=k+0.005
plt.plot(kk,ss,'r*')
plt.xlabel('$k$')
plt.ylabel('$R$',rotation=90)
plt.savefig('模型二.png',dpi=500)
plt.show()
x= [2.09147087e-08 3.08874349e-01 5.14790372e-01 1.40396702e-01
1.52452147e-02 7.72185738e-03]
風(fēng)險(xiǎn)= 0.00772185738406626
【3】模型三
模型三：風(fēng)險(xiǎn)與收益的加權(quán)線性組合。
（1）matlab
clc,clear
M=10000;prob=optimproblem;
x=optimvar('x',6,1,'LowerBound',0);
r=[0.05,0.28,0.21,0.23,0.25];
p=[0,0.01,0.02,0.045,0.065];
q=[0,0.025,0.015,0.055,0.026]';
%w=0:0.1:1;
w=[0.766,0.767,0.810,0.811,0.824,0.825,0.962,0.963,1.0];
V=[];%風(fēng)險(xiǎn)
Q=[];%收益
X=[];%最優(yōu)解
prob.Constraints.con1=(1+p)*x(1:end-1)==M;
prob.Constraints.con2=q(2:end).*x(2:end-1)<=x(end);
for i=1:length(w)
prob.Objective=w(i)*x(end)-(1-w(i))*(r-p)*x(1:end-1);
[sol,fval,flag,out]=solve(prob);
xx=sol.x;
V=[V,max(q.*xx(1:end-1))];
Q=[Q,(r-p)*xx(1:end-1)];
X=[X;xx'];
plot(V,Q,'*-');grid on
xlabel('風(fēng)險(xiǎn)');ylabel('收益')
? ? X
end
V,Q,format
w=（權(quán)重）
? ? [0.766,? ? 0.767,? ? 0.810,? ? 0.811,? ? 0.824,? ? 0.825,? ? 0.962,? ? 0.963,? ? 1.0];
V =（風(fēng)險(xiǎn)）
247.5248 92.2509 92.2509 78.4929 78.4929 59.3960 59.3960 0 0
Q =（收益）
1.0e+03 *
2.6733 2.1648 2.1648 2.1060 2.1060 2.0162 2.0162 0.5000 0.5000
X =
1.0e+04 *
0 0.9901 0 0 0 0.0248
0 0.3690 0.6150 0 0 0.0092
0 0.3690 0.6150 0 0 0.0092
0 0.3140 0.5233 0.1427 0 0.0078
0 0.3140 0.5233 0.1427 0 0.0078
0 0.2376 0.3960 0.1080 0.2284 0.0059（轉(zhuǎn)折點(diǎn)處的X值）
0 0.2376 0.3960 0.1080 0.2284 0.0059
1.0000 0 0 0 0 0
1.0000 0 0 0 0 0
（2）python
from scipy.optimize import linprog
import matplotlib.pyplot as plt
plt.rc('font',size=16)
A=[[0,0.025,0,0,0,-1],[0,0,0.15,0,0,-1],[0,0,0,0.55,0,-1],[0,0,0,0,0.026,-1]]
b=[0,0,0,0]
Aeq=[[1,1.01,1.02,1.045,1.065,0]]
beq=[1]
bound=((0,None),(0,None),(0,None),(0,None),(0,None),(0,None))
s=0;ss=[];aa=[]
while s<=1:
c=[-(1-s)*0.05,-(1-s)*0.27,-(1-s)*0.19,-(1-s)*0.185,-(1-s)*0.185,s]
res=linprog(c,A,b,Aeq,beq)
ss.append(s);aa.append(-res.fun)
s=s+0.01
plt.plot(ss,aa,'r.')
plt.show()import matplotlib.pyplot as plt
import cvxpy as cp
import numpy as np
plt.rc('font',size=16)
x=cp.Variable(6)
a1=np.array([0.025,0.015,0.055,0.026])
a2=np.array([0.05,0.27,0.19,0.185,0.185])
a3=np.array([1,1.01,1.02,1.045,1.065])
con=[cp.multiply(a1,x[1:5])-x[5]<=0,a3@x[:-1]==1,x>=0]
s=0;ss=[];aa=[]
while s<=1:
obj=cp.Minimize(s*x[-1]-(1-s)*a2@x[:-1])
prob=cp.Problem(obj,con)
prob.solve(solver='GLPK_MI')
ss.append(s);aa.append(-prob.value)
s=s+0.01
plt.plot(ss,aa,'r.')
plt.show()
二.整數(shù)規(guī)劃
0.分類
1.通用代碼
%? ? int_LP_1.m
clc,clear,prob=optimproblem;
x=optimvar('x',6,'Type','integer','LowerBound',0);
prob.Objective=sum(x);
% 創(chuàng)建一個(gè)由空優(yōu)化約束組成的 6×1 數(shù)組，使用 constr 初始化用于創(chuàng)建約束表達(dá)式的循環(huán)
coon=optimconstr(6);
a=[35,40,50,45,55,30];
con(1)=x(1)+x(6)>=35;
for i=1:5
con(i+1)=x(i)+x(i+1)>=a(i+1);
end
prob.Constraints.con=con;
[sol,fval,flag]=solve(prob)
x=sol.x
x =
35
5
45
0
55
0intlinprog函數(shù)
x = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)：函數(shù)相比linprog多了一個(gè)參數(shù)intcon，用于標(biāo)定整數(shù)變量的位置：x1、x2為整數(shù)，即intcon = [1 2]
%? ? int_LP_intlinprog.m
c = [-40;-90];
A = [9,7;7,20];
b = [56;70];
[x,fval] = intlinprog(c,[1 2],A,b,[],[],zeros(2,1));
x,maxz=-fval
x =
4.0000
2.0000
maxz =
3402.分支定界法（剪枝）
嘗試一維單調(diào)
%? ? int_branchbunch.m
clear global;
clear;
clc;
global result; % 存儲(chǔ)所有整數(shù)解
global lowerBound; % 下界
global upperBound; % 上界
global count; % 用以判斷是否為第一次分支
count = 1;
f = [-40, -90];
A = [8, 7;7, 20;];
b = [56; 70];
Aeq = [];
beq = [];
lbnd = [0; 0];
ubnd = [inf; inf];
% f = [-3 2 -5];
% A = [1 2 -1;1 4 1;1 1 0;0 4 1];
% b = [2;4;3;6];
% Aeq = [];
% beq = [];
% lbnd=[0 0 0];
% ubnd=[1 1 1];
BinTree = createBinTreeNode({f, A, b, Aeq, beq, lbnd, ubnd});
if ~isempty(result)
[fval,flag]=min(result(:,end)); % result中每一行對(duì)應(yīng)一個(gè)整數(shù)解及對(duì)應(yīng)的函數(shù)值
Result=result(flag,:);
disp('該整數(shù)規(guī)劃問題的最優(yōu)解為：');
disp(Result(1,1:end-1));
disp('該整數(shù)規(guī)劃問題的最優(yōu)值為：');
disp(Result(1,end));
else
disp('該整數(shù)規(guī)劃問題無可行解');
end
function BinTree = createBinTreeNode(binTreeNode)
global result;
global lowerBound;
global upperBound;
global count;
if isempty(binTreeNode)
return;
else
BinTree{1} = binTreeNode;
BinTree{2} = [];
BinTree{3} = [];
[x, fval, exitflag] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, binTreeNode{7});
if count == 1
% upperBound = 0; % 初始下界為空
lowerBound = fval;
count = 2;
end
if ~IsInRange(fval)
return;
end
if exitflag == 1
flag = [];
% 尋找非整數(shù)解分量
for i = 1 : length(x)
if round(x(i)) ~= x(i)
flag = i;
break;
end
end
% 分支
if ~isempty(flag)
lowerBound = max([lowerBound; fval]);
OutputLowerAndUpperBounds();
lbnd = binTreeNode{6};
ubnd = binTreeNode{7};
lbnd(flag, 1) = ceil(x(flag, 1)); % 朝正無窮四舍五入
ubnd(flag, 1) = floor(x(flag, 1));
% 進(jìn)行比較，優(yōu)先選擇目標(biāo)函數(shù)較大的進(jìn)行分支
[~, fval1] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd);
[~, fval2] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7});
if fval1 < fval2
% 創(chuàng)建左子樹
BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});
% 創(chuàng)建右子樹
BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
else
% 創(chuàng)建右子樹
BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
% 創(chuàng)建左子樹
BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});
end
else
upperBound = min([upperBound; fval]);
OutputLowerAndUpperBounds();
result = [result; [x', fval]];
return;
end
else
% 剪枝
return;
end
end
end
function y = IsInRange(fval)
global lowerBound;
global upperBound;
if isempty(upperBound) & fval >= lowerBound
y = 1;
else if (fval >= lowerBound & fval <= upperBound)
y = 1;
else
y = 0;
end
end
end
function y = OutputLowerAndUpperBounds()
global lowerBound;
global upperBound;
disp("此時(shí)下界、上界分別為");
disp(lowerBound);
if isempty(upperBound)
disp(' 無上界')
else
disp(upperBound);
end
end
該整數(shù)規(guī)劃問題的最優(yōu)解為：
4 2
該整數(shù)規(guī)劃問題的最優(yōu)值為：
-3403.割平面法
割去小數(shù)部分，嘗試二維單調(diào)
function [intx,intf] = DividePlane(A,c,b,baseVector)
%功能：用割平面法求解整數(shù)規(guī)劃
%調(diào)用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
%其中, A：約束矩陣；
% c：目標(biāo)函數(shù)系數(shù)向量；
% b：約束右端向量；
% baseVector：初始基向量；
% intx：目標(biāo)函數(shù)取最值時(shí)的自變量值；
% intf：目標(biāo)函數(shù)的最值；
sz = size(A);
nVia = sz(2);%獲取有多少?zèng)Q策變量
n = sz(1);%獲取有多少約束條件
xx = 1:nVia;
if length(baseVector) ~= n
disp('基變量的個(gè)數(shù)要與約束矩陣的行數(shù)相等！');
mx = NaN;
mf = NaN;
return;
end
M = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;
%首先用單純形法求出最優(yōu)解
while 1
[maxs,ind] = max(sigma);
%--------------------用單純形法求最優(yōu)解--------------------------------------
if maxs <= 0 %當(dāng)檢驗(yàn)數(shù)均小于0時(shí)，求得最優(yōu)解。
vr = find(c~=0 ,1,'last');
for l=1:vr
ele = find(baseVector == l,1);
if(isempty(ele))
mx(l) = 0;
else
mx(l)=xb(ele);
end
end
if max(abs(round(mx) - mx))<1.0e-7 %判斷最優(yōu)解是否為整數(shù)解，如果是整數(shù)解。
intx = mx;
intf = mx*c;
return;
else %如果最優(yōu)解不是整數(shù)解時(shí)，構(gòu)建切割方程
sz = size(A);
sr = sz(1);
sc = sz(2);
[max_x, index_x] = max(abs(round(mx) - mx));
[isB, num] = find(index_x == baseVector);
fi = xb(num) - floor(xb(num));
for i=1:(index_x-1)
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
for i=(index_x+1):sc
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
Atmp(1,index_x) = 0; %構(gòu)建對(duì)偶單純形法的初始表格
A = [A zeros(sr,1);-Atmp(1,:) 1];
xb = [xb;-fi];
baseVector = [baseVector sc+1];
sigma = [sigma 0];
%-------------------對(duì)偶單純形法的迭代過程----------------------
while 1
%----------------------------------------------------------
if xb >= 0 %判斷如果右端向量均大于0，求得最優(yōu)解
if max(abs(round(xb) - xb))<1.0e-7
%如果用對(duì)偶單純形法求得了整數(shù)解，則返回最優(yōu)整數(shù)解
vr = find(c~=0 ,1,'last');
for l=1:vr
ele = find(baseVector == l,1);
if(isempty(ele))
mx_1(l) = 0;
else
mx_1(l)=xb(ele);
end
end
intx = mx_1;
intf = mx_1*c;
return;
else %如果對(duì)偶單純形法求得的最優(yōu)解不是整數(shù)解，繼續(xù)添加切割方程
sz = size(A);
sr = sz(1);
sc = sz(2);
[max_x, index_x] = max(abs(round(mx_1) - mx_1));
[isB, num] = find(index_x == baseVector);
fi = xb(num) - floor(xb(num));
for i=1:(index_x-1)
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
for i=(index_x+1):sc
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
Atmp(1,index_x) = 0; %下一次對(duì)偶單純形迭代的初始表格
A = [A zeros(sr,1);-Atmp(1,:) 1];
xb = [xb;-fi];
baseVector = [baseVector sc+1];
sigma = [sigma 0];
continue;
end
else %如果右端向量不全大于0，則進(jìn)行對(duì)偶單純形法的換基變量過程
minb_1 = inf;
chagB_1 = inf;
sA = size(A);
[br,idb] = min(xb);
for j=1:sA(2)
if A(idb,j)<0
bm = sigma(j)/A(idb,j);
if bm<minb_1
minb_1 = bm;
chagB_1 = j;
end
end
end
sigma = sigma -A(idb,:)*minb_1;
xb(idb) = xb(idb)/A(idb,chagB_1);
A(idb,:) = A(idb,:)/A(idb,chagB_1);
for i =1:sA(1)
if i ~= idb
xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
end
end
baseVector(idb) = chagB_1;
end
%------------------------------------------------------------
end
%--------------------對(duì)偶單純形法的迭代過程---------------------
end
else %如果檢驗(yàn)數(shù)有不小于0的，則進(jìn)行單純形算法的迭代過程
minb = inf;
chagB = inf;
for j=1:n
if A(j,ind)>0
bz = xb(j)/A(j,ind);
if bz<minb
minb = bz;
chagB = j;
end
end
end
sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
xb(chagB) = xb(chagB)/A(chagB,ind);
A(chagB,:) = A(chagB,:)/A(chagB,ind);
for i =1:n
if i ~= chagB
xb(i) = xb(i)-A(i,ind)*xb(chagB);
A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
end
end
baseVector(chagB) = ind;
end
M = M + 1;
if (M == 1000000)
disp('找不到最優(yōu)解！');
mx = NaN;
minf = NaN;
return;
end
end
%? ? int_LP_divideplane.m
c = [-1;-1]; % 不要加松弛變量
A = [-1 1 1 0;3 1 0 1]; % 加上松弛變量
b = [1;4];
[x,fval] = DividePlane(A,c,b,[3 4]); % 松弛變量 3 4
x,maxz=fval
x =
1.0000 1.0000
maxz =
-24.匈牙利算法
【1】
【2】
% int_XiongYaLi.m
c = [6 7 11 2;4 5 9 8;3 1 10 4;5 9 8 2];
a = zeros(8,16);
for i = 1:4
a(i,(i-1)*4+1:4*i) = 1;
a(4+i,i:4:16) = 1;
end
b = ones(8,1);
[x,y] = linprog(c,[],[],a,b,zeros(16,1),ones(16,1));
X = reshape(x,4,4)
opt = y
% X =
%
% 0 0 0 1
% 1 0 0 0
% 0 1 0 0
% 0 0 1 0
%
%
% opt =
%
% 15
A=[6 7 11 2;
4 5 9 8;
3 1 10 4;
5 9 8 2];
% A = [3 8 2 10 3;8 7 2 9 7;6 4 2 7 5;8 4 2 3 5;9 10 6 9 10];
B=Hungary(A);
[~,b]=linear_assignment(A,B)
% b =
% 4 1 8 2
function res=Hungary(N)
%輸入的矩陣應(yīng)N*N的
[a,~]=size(N);
%第一步每一行減去當(dāng)前行最小值
for ii = 1:a
N(ii,:)= N(ii,:)-min( N(ii,:));
end
%第二步每一列減去當(dāng)前列最小值
for ii = 1:a
N(:,ii)= N(:,ii)-min( N(:,ii));
end
num=0;
while num~=a
[num,N_min,del_hang,del_lie]=line_count(N);
if num ~=a
for ii=1:a
if del_hang(ii)~=ii
N(ii,:) = N(ii,:)-N_min;
end
if del_lie(ii)==ii
N(:,ii) = N(:,ii)+N_min;
end
end
else
res=N;
end
end
function [num,M_min,del_hang,del_lie]=line_count(M)
[a,~]=size(M);
num=0;
h=0;
del_hang=zeros(a,1);
del_lie=zeros(a,1);
for ii=1:a
del=ii-h;
[~,b]=size(find(M(del,:)==0));
if b>= 2
M(del,:)=[];
h=h+1;
del_hang(ii)=ii; %得到被覆蓋的行數(shù)
num=num+1;
end
end
l=0;
for ii=1:a
del=ii-l;
[b,~]=size(find(M(:,del)==0));
if b >=1
M(:,del)=[];
l=l+1;
del_lie(ii)=ii; %得到被覆蓋的列數(shù)
num=num+1;
end
end
M_min=min(min(M));
end
end
function [place,res]=linear_assignment(M,N)
%N是n維矩陣,N是經(jīng)過Hungary處理的
%M是未處理前的
[a,~]=size(N);
x=0;
place=zeros(1,a);
res=zeros(1,a);
judge=zeros(1,a);
while find(N==0)
for ii=1:a
judge(ii)=length(find(N(ii,:)==0));
end
judge(find(judge==0))=[];
if min(judge)==1
for ii=1:a
if length(find(N(ii,:)==0))==1 %先選出行中只有1個(gè)0
x=x+1;
place(x)=ii+(find(N(ii,:)==0)-1)*a; %得到矩陣中的位置
h=find(N(ii,:)==0);
N(ii,:)=1./zeros(1,a);
N(:,h)=1./zeros(a,1);
end
end
end
for ii=1:a
judge(ii)=length(find(N(ii,:)==0));
end
judge(find(judge==0))=[];
if min(judge)==2
x=x+1;
q=find(N==0);
place(x)=q(1);
N(mod(q(1),a),:)=1./zeros(1,a);
N(:,fix(q(1)/a)+1)=1./zeros(a,1);
end
end
[place,~]=sort(place);
for ii=1:length(place)
res(ii)=M(place(ii));
end
end 5.0-1規(guī)劃
【1】intlinprog函數(shù)
0-1規(guī)劃問題也可以看做區(qū)間在[0,1]的整數(shù)規(guī)劃，下面利用intlinprog函數(shù)進(jìn)行計(jì)算
%? ? intlinprog_01.m
c = [-3 2 -5]';
A = [1 2 -1;1 4 1;1 1 0;0 4 1];
b = [2;4;3;6];
[x,fval] = intlinprog(c,[1 2 3],A,b,[],[],zeros(3,1),ones(3,1));
x,maxz = -fval
x =
1
0
1
maxz =
8【2】指派問題
（1）純0-1整數(shù)規(guī)劃
c=np.array([[4,8,7,15,12],
[7,9,17,14,10],
[6,9,12,8,7],
[6,7,14,6,10],
[6,9,12,10,6]])
x=cp.Variable((5,5),integer=True)
obj=cp.Minimize(cp.sum(cp.multiply(c,x)))
con=[0<=x,x<=1,
? ? cp.sum(x,axis=0,keepdims=True)==1,
cp.sum(x,axis=1,keepdims=True)==1]
prob=cp.Problem(obj,con)
prob.solve(solver='GLPK_MI')
print("最優(yōu)值為：",prob.value)
print("最優(yōu)解為：\n",x.value)
最優(yōu)值為： 34.0
最優(yōu)解為：
[[0. 0. 1. 0. 0.]
[0. 1. 0. 0. 0.]
[1. 0. 0. 0. 0.]
[0. 0. 0. 1. 0.]
[0. 0. 0. 0. 1.]]（2）廣義指派
5*x(1)+4*x(2)<=24+y(1)M;? ? ? ? %M是充分大的數(shù)
7*x(1)+3*x(2)<=45+y(2)M;
y(1)+y(2)=1;i=1,2,3為三種生產(chǎn)方式
x(i)為產(chǎn)量；c(i)為每件產(chǎn)品的成本；k(i)為每種方式的固定成本
min z=(k(1)*y(1)+c(1)*x(1))+(k(2)*y(2)+c(2)*x(2))+(k(3)*y(3)+c(3)*x(3));
y(i)*m<=x(i)<=y(i)*M; %m充分小，M充分大
y(i)=0 or 1;【3】旅行商問題(TSP)
（司P31）比賽項(xiàng)目排序問題中，由于開始項(xiàng)目和結(jié)束項(xiàng)目沒有連接，可引入虛擬項(xiàng)目15，與各項(xiàng)目的權(quán)重為0
三.非線性規(guī)劃
-
算法
【1】二次規(guī)劃
目標(biāo)函數(shù)二次，約束條件線性
H為實(shí)對(duì)稱矩陣
[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)
% DLP_1.m
H=[4,-4;-4,8];
f=[-6;3];
a=[1,1;4,1];
b=[3;9];
[x,fval]=quadprog(H,f,a,b,[],[],zeros(2,1))
x =
2.0854
0.6585
fval =
-5.5976% DLP_2.m
x0 = [1;1];
A=[1,1;4,1];
b=[3;9];
Aeq = [];
beq = [];
lb = [0;0];
[x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb)
function f = fun(x)
f = 2*x(1)-4*x(1)*x(2)+4*x(2)^2-6*x(1)-3*x(2)
end
x =
2.0000
1.0000
fval =
-15.0000【2】無約束非線性規(guī)劃
無約束優(yōu)化問題有局部最優(yōu)解的充分條件：梯度=0；Hesse矩陣正定
%? ? NLP_nocon.m
clc,clear
f = @(x) x(1)^3-x(2)^3+3*x(1)^2+3*x(2)^2-9*x(1);
g = @(x) -f(x);
[m1,n1] = fminunc(f,[0,0])%求極小值
[m2,n2] = fminsearch(g,[0,0]);%求極大值
m2,n2=-n2
m1 =
1.0000 -0.0000
n1 =
-5
m2 =
-3.0000 2.0000
n2 =
31.0000【3】有約束非線性規(guī)劃
%? ? NLP_con_1.m
%[x,fval]=fmincon('fun1',x0,A,b,Aeq,beq,lb,ub,'fun2')
[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')
function f = fun1(x);
f=sum(x.^2)+8;
end
function [g,h] = fun2(x)
g=[-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20];
h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3];
end
x =
0.5522
1.2033
0.9478
y =
10.6511% NLP_con_2.m
clc,clear
x = optimvar('x',3,'LowerBound',0);
prob = optimproblem('Objective',sum(x.^2)+8);
con1 = [-x(1)^2+x(2)-x(3)^2 <= 0
x(1)+x(2)^2+x(3)^3 <= 20];
con2 = [-x(1)-x(2)^2+2 == 0
x(2)+2*x(3)^2 == 3];
prob.Constraints.con1 = con1;
prob.Constraints.con2 = con2;
x0.x = rand(3,1);
[s,f,flag,out] = solve(prob,x0);
s.x,f（1）只有等式約束
拉格朗日法
（2）一般形式
罰函數(shù)法
不等式約束g(x,i)<=0? ? 等式約束h(x,j)=0
不等式 g(x,i)<=0 轉(zhuǎn)化為 等式 max(0,g(x,i))=0
罰函數(shù)：T(x.M)=f(x)+M*sum(max(0，g(x,i)))+M*sum(h(x,j)^2)【4】凸規(guī)劃（f''>=0）
（1）判定：Hesse處處半正定
%? ? NLP_TU.m
clc,clear
prob=optimproblem;
x=optimvar('x',2,'LowerBound',0);
prob.Objective=sum(x.^2)-4*x(1)+4;
con=[-x(1)+x(2)-2<=0
x(1)^2-x(2)+1<=0];
prob.Constraints.con=con;
x0.x=rand(2,1)
[sol,fval,flag,out]=solve(prob,x0),sol.x（2）K-T條件
只要是最優(yōu)點(diǎn)，一定滿足K-T條件
若為凸規(guī)劃，則充要條件
【5】蒙特卡洛法
%? ? MTKL.m
clc,clear
rng(0) %rng('shuffle')
p0=0;n=10^6;tic
for i=1:n
x=randi([0,99],1,5);
[f,g]=mengte(x);
if all(g<=0)
if p0<f
x0=x;p0=f;
end
end
end
x0,p0,toc
function [f,g]=mengte(x);
f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)^2-8*x(1)-2*x(2)-3*x(3)-x(4)-2*x(5);
g=[sum(x)-400
x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800
2*x(1)+x(2)+6*x(3)-200
x(3)+x(4)+5*x(5)-200];
end
x0 =
46 98 1 99 3
p0 =
502732.LINGO
model:
sets:
row/1..4/:b;
col/1..5/:c1,c2,x;
link(row,col):a;
endsets
data:
c1=1,1,3,4,2;
c2=-8,-2,-3,-1,-2;
a=1 1 1 1 1
1 2 2 1 6
2 1 6 0 0
0 0 1 1 5;
b=400,800,200,200;
enddata
max=@sum(col:c1*x^2+c2*x);
@for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));
@for(col:@gin(x));
@for(col:@bnd(0,x,99));
end
3.案例+靈敏度分析
（司P42）
clc,clear,format long g,close all
syms x1 x2
f=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2);
f=simplify(f);
f1=diff(f,x1);f2=diff(f,x2);
[x10,x20]=solve(f1,f2);
x10=round(double(x10));x20=round(double(x20));
f0=subs(f,{x1,x2},{x10,x20});
f0=double(f0);
subplot(121),fmesh(f,[0,10000,0,10000]),title('')
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex')
subplot(122),L=fcontour(f,[0,10000,0,10000]);
contour(L.XData,L.YData,L.ZData,'ShowText','on')
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
p1=339-0.01*x10-0.003*x20
p2=399-0.004*x10-0.01*x20
c=400000+195*x10+225*x20 %總支出
rate=f0/c %利潤
x10
x20
靈敏度分析
clc,clear,format long g,close all
syms x1 x2 a
f=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2);
f=simplify(f);
f1=diff(f,x1);f2=diff(f,x2);
[x10,x20]=solve(f1,f2);
subplot(121),fplot(x10,[0.002,0.02]),title('') %x1關(guān)于a的曲線
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
subplot(122),fplot(x20,[0.002,0.02]),title('') %x2關(guān)于a的曲線
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
dx1=diff(x10,a);dx10=subs(dx1,a,0.01);dx10=double(dx10);
sx1a=dx10*0.01/4735;
dx2=diff(x20,a);dx20=subs(dx2,a,0.01);dx20=double(dx20);
sx2a=dx20*0.01/7043;
F=subs(f,{x1,x2},{x10,x20});
F=simplify(F);
figure,fplot(F,[0.002,0.02]),title('')
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
Sya=-4735^2*0.01/553641;
f3=subs(f,{x1,x2,a},{4735,7043,0.011});
f3=double(f3);
f4=subs(F,a,0.011) %近似最優(yōu)利潤
f4=double(f4) %最優(yōu)利潤
delta=(f4-f3)/f4 %利潤的相對(duì)誤差四.多目標(biāo)規(guī)劃
1.多目標(biāo)規(guī)劃
max f1(x)=2x(1)+3x(2)
min f2(x)=x(1)+2x(2)
0.5x(1)+0.25x(2)<=8
0.2x(1)+0.2x(2)<=4
x(1)+x(2)>=10
x(1)>=0
x(2)>=0
【1】線性加權(quán)
min 0.5(-f1)+0.5(f2)
結(jié)果：
x(1)=7，x(2)=13
【2】理想點(diǎn)
分別求解得f1=-53，f2=10
新的目標(biāo)函數(shù) min f=[-2x(1)-3x(2)+53]^2+[x(1)+2x(2)-10]^2
結(jié)果：
x(1)=13.36，x(2)=5.28
【3】序貫
分別求解得f1=-53
加入條件-2x(1)-3x(2)=-53求f2
結(jié)果：
x(1)=7，x(2)=13
clc,clear,prob=optimproblem;
x=optimvar('x',2,'LowerBound',0);
c1=[-2,-3];c2=[1,2];
a=[0.5,0.25;0.2,0.2;1,5;-1,-1];
b=[8;4;72;-10];
prob.Constraints.con1=a*x<=b
obj1=0.5*c1*x+0.5*c2*x
prob1=prob;prob1.Objective=obj1;
[sol1,fval1]=solve(prob1),sx=sol1.x
f1=-c1*sx,f2=c2*sx
prob21=prob;prob21.Objective=c1*x;
[sol21,fval21]=solve(prob21),sx21=sol21.x
prob22=prob;prob22.Objective=c2*x;
[sol22,fval22]=solve(prob22),sx22=sol22.x
prob23=prob;
prob23.Objective=(c1*x-fval21)^2+(c2*x-fval22)^2;
[sol23,fval23]=solve(prob23),sx23=sol23.x
prob3=prob;prob3.Objective=c2*x;
prob3.Constraints.con2=c1*x==fval21;
[sol3,fval3]=solve(prob3),sx3=sol3.x2.目標(biāo)規(guī)劃
正負(fù)偏差變量
軟硬約束
“盡量”
優(yōu)先權(quán)
序貫算法：文章來源：http://www.zghlxwxcb.cn/news/detail-649112.html
根據(jù)優(yōu)先級(jí)排序，將目標(biāo)規(guī)劃分解成一系列單目標(biāo)規(guī)劃，依次求解文章來源地址http://www.zghlxwxcb.cn/news/detail-649112.html
clc,clear
x=optimvar('x',2,'LowerBound',0);
dp=optimvar('dp',4,'LowerBound',0);
dm=optimvar('dm',4,'LowerBound',0);
p=optimproblem('ObjectiveSense','min');
p.Constraints.con1=2*sum(x)<=12;
con2=[200*x(1)+300*x(2)+dm(1)-dp(1)==1500
2*x(1)-x(2)+dm(2)-dp(2)==0
4*x(1)+dm(3)-dp(3)==16
5*x(2)+dm(4)-dp(4)==15];
p.Constraints.con2=con2;
goal=100000*ones(3,1);
mobj=[dm(1);dp(2)+dm(2);3*dp(3)+3*dm(3)+dp(4)];
for i=1:3
p.Constraints.con3=[mobj<=goal];
p.Objective=mobj(i);
[sx,fval]=solve(p);
fprintf('第%d級(jí)目標(biāo)計(jì)算結(jié)果如下:\n',i)
fval,xx=sx.x,sdm=sx.dm,sdp=sx.dp
goal(i)=fval;
end
到了這里，關(guān)于數(shù)學(xué)建?！€性規(guī)劃類的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
