matlab进阶

大家好，欢迎来到IT知识分享网。

设置matlab路径

你写了一个外部函数，想在其他目录下直接拿来用的话需要设置下工作路径

 把一些含有需要的外部函数文件的文件夹添加就好了

 线性

ＬＰ　ｌｉｎｅａｅｒ　ｐｒｏｇｒａｍｉｎｇ线性规划

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．ｍ和intprog．ｍ配置到路径

branchbound．ｍ

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．ｍ

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.ｍ　将这个.ｍ文件配置到路径里面

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 == ) 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])

别忘了ｃ里是负的。

0,1选择问题，互斥约束问题

生产线生产产品，新旧工艺只能用一个，用了新工艺就不能用旧工艺；用了旧工艺就不能用新工艺。

去很大很大的值。当时，执行第一个方程,因为是很大的数，让第二个方程没有什么意义了。

都是无穷大的数