设置matlab路径
你写了一个外部函数,想在其他目录下直接拿来用的话需要设置下工作路径
把一些含有需要的外部函数文件的文件夹添加就好了
线性
LP lineaer programing线性规划
matlab只能求最小,求最大需要加负号
Aeq 就是A equal
当线性规划的约束中有aiXi=bi的条件时,用Aeq和beq来保证等号的成立,就是说aiXi=bi那个xi对应的Aeq位子取ai,其他取0,beq取值是有几个aiXi=bi就取几个bi。
例子
matlab只能求最小,求最大需要加负号
除了等式不需要加负号,其他都要加。
这里的f是目标函数的,a.b 是不等式的。aeq,beq是等式的。用linprog求解,x,y=-y
求不出来,因为没具体的数。看思想。
多目标规划改成线性规划。
投资问题
是风险损失率,
是风险额,总体风险率最小。
是最大金额
总风险减去总收益最小就可以。或者是总收益减去总风险最大就可以。
假设以a 方案构建
这些数字是用已知的图表得到的。其中
是5%
clc ,clear
a=0,hold on
while a<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=a*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(a,Q,'*k'); %k是黑色颜色的代号
a=a+0.001;
end
xlabel('a'),ylabel('Q')
函数图形用黑色颜色和“*”符号来表示a和Q的函数曲线关系。
这里,plot是绘图函数,a和Q是一组数据,k是黑色颜色的代号,“*”在图形中表示a和Q的对应值。
整数规划(Integer Programming,IP)
(94条消息) 整数规划之分支定界法_咖瑞芝的博客-CSDN博客_整数规划分支定界法
(94条消息) 整数规划:分支定界法_胡拉哥-CSDN博客_整数规划的分支定界法
分支定界法
分支定界法的核心思想就是分枝和剪枝。当我们不考虑所求解必须是整数这个条件时,用单纯形法可求出最优解,但是这个解往往不全是整数,因此我们采用剪枝的方式一点一点缩小范围,直到所求解为整数解。
剪枝:如果某一个子问题无可行解或者最优值小于原来的下界,则称这个分支已经查清,将该支剪掉,不再计算。
解出现小数的形式,取上界和下界。也就是说,假如x1=2.25和x2=3.75.先选一个分支x2,加新的约束条件。这是因为不确定整数解到底在那个方向能取到最优值,因此需要考虑两种情况。
求解整数规划的分支定界算法是一个树搜索(Tree-Search)算法。它的基本思想是考虑整数规划问题的松弛问题,通过增加“整性”约束条件来构造子问题,求解所有子问题从而得到最优的整数解。通过“定界”,可以避免不必要的搜索。
将branchbound.m和intprog.m配置到路径
branchbound.m
function [newx,newfval,status,newbound] = branchbound(f,A,B,I,x,fval,bound,Aeq,Beq,lb,ub,e)
% 分支定界法求解整数规划
% f,A,B,Aeq,Beq,lb,ub与线性规划相同
% I为整数限制变量的向量
% x为初始解,fval为初始值
options = optimset('display','off');
[x0,fval0,status0]=linprog(f,A,B,Aeq,Beq,lb,ub,[],options);
%递归中的最终退出条件
%无解或者解比现有上界大则返回原解
if status0 <= 0 || fval0 >= bound
newx = x;
newfval = fval;
newbound = bound;
status = status0;
return;
end
%是否为整数解,如果是整数解则返回
intindex = find(abs(x0(I) - round(x0(I))) > e);
if isempty(intindex) %判断是否为空值
newx(I) = round(x0(I));
newfval = fval0;
newbound = fval0;
status = 1;
return;
end
%当有非整可行解时,则进行分支求解
%此时必定会有整数解或空解
%找到第一个不满足整数要求的变量
n = I(intindex(1));
addA = zeros(1,length(f));
addA(n) = 1;
%构造第一个分支 x<=floor(x(n))
A = [A;addA];
B = [B,floor(x(n))];%向下取整
[x1,fval1,status1,bound1] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
A(end,:) = [];
B(:,end) = [];
%解得第一个分支,若为更优解则替换,若不是则保持原状
status = status1;
if status1 > 0 && bound1 < bound
newx = x1;
newfval = fval1;
bound = fval1;
newbound = bound1;
else
newx = x0;
newfval = fval0;
newbound = bound;
end
%构造第二分支
A = [A;-addA];
B = [B,-ceil(x(n))];%向上取整
[x2,fval2,status2,bound2] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
A(end,:) = [];
B(:,end) = [];
%解得第二分支,并与第一分支做比较,如果更优则替换
if status2 > 0 && bound2 < bound
status = status2;
newx = x2;
newfval = fval2;
newbound = bound2;
end
intprog.m
function [x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
%整数规划求解函数 intprog()
% 其中 f为目标函数向量
% A和B为不等式约束 Aeq与Beq为等式约束
% I为整数约束
% lb与ub分别为变量下界与上界
% x为最优解,fval为最优值
%例子:
% maximize 20 x1 + 10 x2
% S.T.
% 5 x1 + 4 x2 <=24
% 2 x1 + 5 x2 <=13
% x1, x2 >=0
% x1, x2是整数
% f=[-20, -10];
% A=[ 5 4; 2 5];
% B=[24; 13];
% lb=[0 0];
% ub=[inf inf];
% I=[1,2];
% e=0.000001;
% [x v s]= IP(f,A,B,I,[],[],lb,ub,,e)
% x = 4 1 v = -90.0000 s = 1
% 控制输入参数
if nargin < 9, e = 0.00001;
if nargin < 8, ub = [];
if nargin < 7, lb = [];
if nargin < 6, Beq = [];
if nargin < 5, Aeq = [];
if nargin < 4, I = [1:length(f)];
end, end, end, end, end, end
%求解整数规划对应的线性规划,判断是否有解
options = optimset('display','off');
[x0,fval0,exitflag] = linprog(f,A,B,Aeq,Beq,lb,ub,[],options);
if exitflag < 0
disp('没有合适整数解');
x = x0;
fval = fval0;
status = exitflag;
return;
else
%采用分支定界法求解
bound = inf;
[x,fval,status] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
end
%例子1
% f = [-40 -90];
%A = [9 7;7 20];
%B = [56 70];
% lb = [0 0]';
%例子2
f = [-20 -10];
A = [5 4;2 5];
B = [24 13];
lb = [0 0];
[x,fval,status] = intprog(f,A,B,[1 2],[],[],lb)
branchbound.m 和intprog.m都是拿来直接用的。用来求解整数问题。
割平面算法
首先不考虑变量 xi是整数这一条件,但增加线性约束条件( 用几何术语,称为割平面) 使得由原可行域中切割掉一部分,这部分只包含非整数解, 但没有切割掉任何整数可行解。这个方法就是指出怎样找到适当的割平面(不见得一次就找到) ,使切割后最终得到这样的可行域, 它的一个有整数坐标的极点恰好是问题的最优解。
有几个不等式就要引入几个松弛变量。
DividePlane.m 将这个.m文件配置到路径里面
function [intx,intf] = DividePlane(A,c,b,baseVector)
%功能:用割平面法求解整数规划
%调用格式:[intx,intf]=DividePlane(A,c,b,baseVector)
%其中,A:约束矩阵;
% c:目标函数系数向量;
% b:约束右端向量;
% baseVector:初始基向量;
% intx:目标函数取最小值时的自变量值;
% intf:目标函数的最小值;
sz = size(A);
nVia = sz(2);%获取有多少决策变量
n = sz(1);%获取有多少约束条件
xx = 1:nVia;
if length(baseVector) ~= n
disp('基变量的个数要与约束矩阵的行数相等!');
mx = NaN;
mf = NaN;
return;
end
M = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;
%首先用单纯形法求出最优解
while 1
[maxs,ind] = max(sigma);
%--------------------用单纯形法求最优解--------------------------------------
if maxs <= 0 %当检验数均小于0时,求得最优解。
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 %判断最优解是否为整数解,如果是整数解。
intx = mx;
intf = mx*c;
return;
else %如果最优解不是整数解时,构建切割方程
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; %构建对偶单纯形法的初始表格
A = [A zeros(sr,1);-Atmp(1,:) 1];
xb = [xb;-fi];
baseVector = [baseVector sc+1];
sigma = [sigma 0];
%-------------------对偶单纯形法的迭代过程----------------------
while 1
%----------------------------------------------------------
if xb >= 0 %判断如果右端向量均大于0,求得最优解
if max(abs(round(xb) - xb))<1.0e-7 %如果用对偶单纯形法求得了整数解,则返回最优整数解
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 %如果对偶单纯形法求得的最优解不是整数解,继续添加切割方程
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; %下一次对偶单纯形迭代的初始表格
A = [A zeros(sr,1);-Atmp(1,:) 1];
xb = [xb;-fi];
baseVector = [baseVector sc+1];
sigma = [sigma 0];
continue;
end
else %如果右端向量不全大于0,则进行对偶单纯形法的换基变量过程
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
%--------------------对偶单纯形法的迭代过程---------------------
end
else %如果检验数有不小于0的,则进行单纯形算法的迭代过程
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('找不到最优解!');
mx = NaN;
minf = NaN;
return;
end
end
A=[0.01 0.01 0.01 0.03 0.03 0.03 1 0 0 0;0.02 0 0 0.05 0 0 0 1 0 0;0 0.02 0 0 0.05 0 0 0 1 0;0 0 0.03 0 0 0.08 0 0 0 1];
c=[-20;-14;-16;-36;-32;-30];
b=[850;700;100;900];
[intx,intf]=DividePlane(A,c,b,[7 8 9 10])
别忘了c里是负的。
0,1选择问题,互斥约束问题
生产线生产产品,新旧工艺只能用一个,用了新工艺就不能用旧工艺;用了旧工艺就不能用新工艺。
去很大很大的值。当
时,执行第一个方程
,因为
是很大的数,让第二个方程没有什么意义了。
都是无穷大的数