如何构造求解规划模型——数学建模（MATLAB）

最新推荐文章于 2024-06-18 08:30:00 发布

Zg·ln

最新推荐文章于 2024-06-18 08:30:00 发布

阅读量966

点赞数 5

分类专栏：数学建模文章标签：数学建模 matlab

本文链接：https://blog.csdn.net/Zg_ln/article/details/133744961

版权

数学建模专栏收录该内容

1 篇文章 0 订阅

订阅专栏

规划模型主要分为线性规划、整数规划、混合整数规划、非线性规划、多目标规划。接下来会一一讲述各类使用方法：可以直接用！！！

线性规划

线性规划就是在一组线性约束条件下,求线性目标函数的最大或最小值。"线性"意味着所有变量都是一次方。

使用场景：题目中提到"怎样安排/分配”"尽量多(少)”“最多(少)” "利润最大”“最合理”等词;

%例子:

% max 2 x1 + 3 x2 -5x3

% S.T.

% -2x1 + 5 x2 -x3<=-10

% x1 + 3 x2 +x3<=12

% x1 + x2 + x3=7

% x1, x2 ,x3>=0

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;

在MATLAB中解决线性规划问题，一般使用linprog函数，其中f为目标函数，a,b为约束条件，aeq,beq为决策变量的范围即上下界，zeros(3,1)为决策变量初始值

整数规划

整数规划所求解是整数与线性规划所求的实数不同，依据不同的问题会有不同的方法求解，主要有分支定界法、割平面法、隐枚举法、匈牙利法、蒙特卡洛法。分支定界法、割平面法适用于一般整数规划，隐枚举法、匈牙利法适用于0-1规划问题、指派问题，蒙特卡洛法都适用，但求的解为近似值。

分支定界法

%例子:

% 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]= intprog(f,A,B,I,[],[],lb,ub,,e)

使用intprog函数进行求解，intprog函数并不在MATLAB函数库中，在下面给出：

代码：

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为最优值

% 控制输入参数

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

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

割平面法

%例子：

%max Z =_20x1 +14x2 +16x3 +36x4 + 32x5 +30x6

%0.01x1 +0.01x2 +0.01x3 +0.03x4 +0.03x5+ 0.03x6 <= 850

%0.02x1 +0.05x4 <= 700

%0.02x1+0.05x5<=100

%0.03x3+0.08x6<=900

%Xi>=0,且为整数，i=1,...,6

化不等式为等式

%min -20x1 -14x2 -16x3 -36x4 - 32x5 -30x6

%0.01x1 +0.01x2 +0.01x3 +0.03x4 +0.03x5+ 0.03x6 +x7= 850

%0.02x1 +0.05x4 +x8= 700

%0.02x1+0.05x5 +x9100

%0.03x3+0.08x6 +x10=900

%Xi>=0,且为整数，i=1,...,10

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];

b = [850;700;100;900];

c = [-20;-14;-16;-36;-32;-30];

[intx,intf] = DividePlane(A,c,b,[7 8 9 10])

使用DividePlane函数进行求解，DividePlane函数并不在MATLAB函数库中，在下面给出：

代码：

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

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

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

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

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

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

baseVector(chagB) = ind;

end

M = M + 1;

if (M == 1000000)

disp('找不到最优解！');

mx = NaN;

minf = NaN;

return;

end

隐枚举法

%例子：

%min -3 x1+2 x2-5 x3

% x1+2 x2- x3<=2;

% x1+4 x2+ x3<=4;

% x1+ x2<=3;

%4 x1+x3<=6;

%x1,x2,x3=0,1;

A=[1 2 -1;1 4 1;1 1 0;4 0 1];

b=[2;4;3;6];

c=[-3 2 -5];

x=solve01ILP(c, A, b)

使用solve01ILP函数进行求解，solve01ILP函数并不在MATLAB函数库中，在下面给出：

代码：

solve01ILP函数

function x = solve01ILP(c, A, b)

%在这段代码中，c是目标函数的系数向量，A和b分别是约束矩阵和约束向量。

%函数solve01ILP通过调用intlinprog求解0-1整数线性规划问题。函数返回一个包含0-1解的向量x

n = size(A,2);

lb = zeros(n, 1);

ub = ones(n, 1);

intcon = 1:n;

options = optimoptions('intlinprog', 'Display', 'off');

x = intlinprog(c, intcon, A, b, [], [], lb, ub, options);

end

匈牙利法

例子：

有n项不同的任务，需要n个人分别完成其中的1项，每个人完成任务的时间不一样。于是就有一个问题，如何分配任务使得花费时间最少。

%C为关于任务与人数n*n的矩阵

C=[2 15 13 4;10 4 14 15;9 14 16 13;7 8 11 9];

[Matching,Cost]=Hungarian (C);

disp('最优解矩阵为:');%输出指派方案和最优值

Matching

disp('最优解为:');

Cost

使用Hungarian函数进行求解，Hungarian函数并不在MATLAB函数库中，在下面给出：

注意：若出现m*n的矩阵，可以补全成为方阵（补一个最小值或一个最大值）

代码：

Hungarian函数

%指派问题的匈牙利算法，输入矩阵，a(ij)为i指派给j,第i人干第j个工作

function [Matching,Cost] = Hungarian(Perf)

% [MATCHING,COST] = Hungarian_New(WEIGHTS)

% A function for finding a minimum edge weight matching given a MxN Edge

% weight matrix WEIGHTS using the Hungarian Algorithm.

% An edge weight of Inf indicates that the pair of vertices given by its

% position have no adjacent edge.

% MATCHING return a MxN matrix with ones in the place of the matchings and

% zeros elsewhere.

% COST returns the cost of the minimum matching

% Written by: Alex Melin 30 June 2006

% Initialize Variables

Matching = zeros(size(Perf));

% Condense the Performance Matrix by removing any unconnected vertices to

% increase the speed of the algorithm

% Find the number in each column that are connected

num_y = sum(~isinf(Perf),1);

% Find the number in each row that are connected

num_x = sum(~isinf(Perf),2);

% Find the columns(vertices) and rows(vertices) that are isolated

x_con = find(num_x~=0);

y_con = find(num_y~=0);

% Assemble Condensed Performance Matrix

P_size = max(length(x_con),length(y_con));

P_cond = zeros(P_size);

P_cond(1:length(x_con),1:length(y_con)) = Perf(x_con,y_con);

if isempty(P_cond)

Cost = 0;

return

end

% Ensure that a perfect matching exists

% Calculate a form of the Edge Matrix

Edge = P_cond;

Edge(P_cond~=Inf) = 0;

% Find the deficiency(CNUM) in the Edge Matrix

cnum = min_line_cover(Edge);

% Project additional vertices and edges so that a perfect matching

% exists

Pmax = max(max(P_cond(P_cond~=Inf)));

P_size = length(P_cond)+cnum;

P_cond = ones(P_size)*Pmax;

P_cond(1:length(x_con),1:length(y_con)) = Perf(x_con,y_con);

%*************************************************

% MAIN PROGRAM: CONTROLS WHICH STEP IS EXECUTED

%*************************************************

exit_flag = 1;

stepnum = 1;

while exit_flag

switch stepnum

case 1

[P_cond,stepnum] = step1(P_cond);

case 2

[r_cov,c_cov,M,stepnum] = step2(P_cond);

case 3

[c_cov,stepnum] = step3(M,P_size);

case 4

[M,r_cov,c_cov,Z_r,Z_c,stepnum] = step4(P_cond,r_cov,c_cov,M);

case 5

[M,r_cov,c_cov,stepnum] = step5(M,Z_r,Z_c,r_cov,c_cov);

case 6

[P_cond,stepnum] = step6(P_cond,r_cov,c_cov);

case 7

exit_flag = 0;

end

% Remove all the virtual satellites and targets and uncondense the

% Matching to the size of the original performance matrix.

Matching(x_con,y_con) = M(1:length(x_con),1:length(y_con));

Cost = sum(sum(Perf(Matching==1)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% STEP 1: Find the smallest number of zeros in each row

% and subtract that minimum from its row

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [P_cond,stepnum] = step1(P_cond)

P_size = length(P_cond);

% Loop throught each row

for ii = 1:P_size

rmin = min(P_cond(ii,:));

P_cond(ii,:) = P_cond(ii,:)-rmin;

end

stepnum = 2;

%**************************************************************************

% STEP 2: Find a zero in P_cond. If there are no starred zeros in its

% column or row start the zero. Repeat for each zero

%**************************************************************************

function [r_cov,c_cov,M,stepnum] = step2(P_cond)

% Define variables

P_size = length(P_cond);

r_cov = zeros(P_size,1); % A vector that shows if a row is covered

c_cov = zeros(P_size,1); % A vector that shows if a column is covered

M = zeros(P_size); % A mask that shows if a position is starred or primed

for ii = 1:P_size

for jj = 1:P_size

if P_cond(ii,jj) == 0 && r_cov(ii) == 0 && c_cov(jj) == 0

M(ii,jj) = 1;

r_cov(ii) = 1;

c_cov(jj) = 1;

end

% Re-initialize the cover vectors

r_cov = zeros(P_size,1); % A vector that shows if a row is covered

c_cov = zeros(P_size,1); % A vector that shows if a column is covered

stepnum = 3;

%**************************************************************************

% STEP 3: Cover each column with a starred zero. If all the columns are

% covered then the matching is maximum

%**************************************************************************

function [c_cov,stepnum] = step3(M,P_size)

c_cov = sum(M,1);

if sum(c_cov) == P_size

stepnum = 7;

else

stepnum = 4;

end

%**************************************************************************

% STEP 4: Find a noncovered zero and prime it. If there is no starred

% zero in the row containing this primed zero, Go to Step 5.

% Otherwise, cover this row and uncover the column containing

% the starred zero. Continue in this manner until there are no

% uncovered zeros left. Save the smallest uncovered value and

% Go to Step 6.

%**************************************************************************

function [M,r_cov,c_cov,Z_r,Z_c,stepnum] = step4(P_cond,r_cov,c_cov,M)

P_size = length(P_cond);

zflag = 1;

while zflag

% Find the first uncovered zero

row = 0; col = 0; exit_flag = 1;

ii = 1; jj = 1;

while exit_flag

if P_cond(ii,jj) == 0 && r_cov(ii) == 0 && c_cov(jj) == 0

row = ii;

col = jj;

exit_flag = 0;

end

jj = jj + 1;

if jj > P_size; jj = 1; ii = ii+1; end

if ii > P_size; exit_flag = 0; end

end

% If there are no uncovered zeros go to step 6

if row == 0

stepnum = 6;

zflag = 0;

Z_r = 0;

Z_c = 0;

else

% Prime the uncovered zero

M(row,col) = 2;

% If there is a starred zero in that row

% Cover the row and uncover the column containing the zero

if sum(find(M(row,:)==1)) ~= 0

r_cov(row) = 1;

zcol = find(M(row,:)==1);

c_cov(zcol) = 0;

else

stepnum = 5;

zflag = 0;

Z_r = row;

Z_c = col;

end

%**************************************************************************

% STEP 5: Construct a series of alternating primed and starred zeros as

% follows. Let Z0 represent the uncovered primed zero found in Step 4.

% Let Z1 denote the starred zero in the column of Z0 (if any).

% Let Z2 denote the primed zero in the row of Z1 (there will always

% be one). Continue until the series terminates at a primed zero

% that has no starred zero in its column. Unstar each starred

% zero of the series, star each primed zero of the series, erase

% all primes and uncover every line in the matrix. Return to Step 3.

%**************************************************************************

function [M,r_cov,c_cov,stepnum] = step5(M,Z_r,Z_c,r_cov,c_cov)

zflag = 1;

ii = 1;

while zflag

% Find the index number of the starred zero in the column

rindex = find(M(:,Z_c(ii))==1);

if rindex > 0

% Save the starred zero

ii = ii+1;

% Save the row of the starred zero

Z_r(ii,1) = rindex;

% The column of the starred zero is the same as the column of the

% primed zero

Z_c(ii,1) = Z_c(ii-1);

else

zflag = 0;

end

% Continue if there is a starred zero in the column of the primed zero

if zflag == 1;

% Find the column of the primed zero in the last starred zeros row

cindex = find(M(Z_r(ii),:)==2);

ii = ii+1;

Z_r(ii,1) = Z_r(ii-1);

Z_c(ii,1) = cindex;

end

% UNSTAR all the starred zeros in the path and STAR all primed zeros

for ii = 1:length(Z_r)

if M(Z_r(ii),Z_c(ii)) == 1

M(Z_r(ii),Z_c(ii)) = 0;

else

M(Z_r(ii),Z_c(ii)) = 1;

end

% Clear the covers

r_cov = r_cov.*0;

c_cov = c_cov.*0;

% Remove all the primes

M(M==2) = 0;

stepnum = 3;

% *************************************************************************

% STEP 6: Add the minimum uncovered value to every element of each covered

% row, and subtract it from every element of each uncovered column.

% Return to Step 4 without altering any stars, primes, or covered lines.

%**************************************************************************

function [P_cond,stepnum] = step6(P_cond,r_cov,c_cov)

a = find(r_cov == 0);

b = find(c_cov == 0);

minval = min(min(P_cond(a,b)));

P_cond(find(r_cov == 1),:) = P_cond(find(r_cov == 1),:) + minval;

P_cond(:,find(c_cov == 0)) = P_cond(:,find(c_cov == 0)) - minval;

stepnum = 4;

function cnum = min_line_cover(Edge)

% Step 2

[r_cov,c_cov,M,stepnum] = step2(Edge);

% Step 3

[c_cov,stepnum] = step3(M,length(Edge));

% Step 4

[M,r_cov,c_cov,Z_r,Z_c,stepnum] = step4(Edge,r_cov,c_cov,M);

% Calculate the deficiency

cnum = length(Edge)-sum(r_cov)-sum(c_cov);

蒙特卡洛法

适用于任何情况，但是近似值。其运行过程主要在脚本中

代码：

蒙特卡洛函数

function [f,g] = mengte(x)

% mengte()蒙特卡洛函数输入对应解的向量，返回目标函数值f 以及约束函数值g

% 目标函数

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 % sum()函数用于对向量求和，或对矩阵的列求和，即x(1) + x(2) + x(3) + x(4) + x(5)

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

% 蒙特卡洛法

% 用matlab来进行蒙特卡洛法，可以得到一个满意的解，理论上也是非常接近精确解

% 但它并不是精确解，所以它每一次的解都会不一样

rand('state', sum(clock)); % 初始化随机数发生器

% 用rand('state',X)来设置随机数流的状态，使随机数随X来改变

% 与C++中的seed相似

p0=0;

tic % 计时开始呼应上文的clock

for i = 1:10^6 %随机次数，越大精度越高

x=randi([0,99], 1, 5); % 产生一行五列的在区间[0, 99]上的随机整数

% randi()函数，用于产生一个随机整数矩阵

[f, g] = mengte(x); % 调用mengte()函数

if all(g<=0) % 使用约束函数进行限制（可调）

% all()函数对向量/矩阵所有项进行相同的逻辑比较，并返回bool值

if p0 < f % 用p0储存上一个的f值，并通过比较来得到最大的f

x0 = x; % 记录目前达到最大值的x值

p0 = f; % 记录目前的最大值

end

x0,p0 % 输出x0和p0的值

混合整数规划

使用MATLAB自带的intlinprog函数，注意在这里多了intcon 表示整数所在的位置

代码：

%intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)

% f ：目标函数的系数矩阵

% intcon ：整数所在位置

% A ：不等式约束的变量系数矩阵

% b ：不等式约束的资源数

% Aeq ：等式约束的变量系数矩阵

% beq ：等式约束的资源数

% lb ：变量约束下限

% ub ：变量约束上限

f = [7 5 9 6 3];

ic = [1,2,3,4,5];

A = [56,20,54,42,15;1,4,1,0,0;-1,-2,0,-1,-2];

b = [100;4;-2];

lb = zeros(5,1); % 生成5×1的 0 矩阵

ub = ones(5,1); % 生成5×1的 1 矩阵

[x,fval]=intlinprog(f,ic,A,b,[],[],lb,ub)

非线性规划

一般使用MATLAB自带的fmincon函数

代码：

脚本

[x,y] = fmincon('fun1',[0;0],[],[],[],[],[0;0;0],[],'fun2')

fun1函数

function f = fun1(x)

f = sum(x.^2) + 8;（可修改）

%目标函数

fun2函数

function [g,h]=fun2(x)

g = [-x(1)^2+x(2)-x(3)^2（可修改）

x(1)+x(2)^2+x(3)^3-20] %g非线性不等式

h=[-x(1)-x(2)^2+2 （可修改）

x(2)+2*x(3)^2-3] %h是非线性等式

多目标规划

多目标规划转化为单日标规划划方法:

评价函数法

功效系数法

约束法

分层序列法

小极大法

多目标规划的MATLAB 函数 fgoalattain 的用法为

[x,fval]=fgoalattain (“fun',x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)

例子：

脚本（可修改）

a=[-1 -1 0 0;0 0 -1 -1;3 0 2 0;0 3 0 2]

b=[-30,-30,120,48]

c1=-[100,90,80,70]

c2=[0,3,0,2]

[x1,goal1]=linprog(c1,a,b,[],[],zeros(4,1))

[x2,goal2]=linprog(c2,a,b,[],[],zeros(4,1))

goal=[goal1;goal2]

[x,fval]=fgoalattain('Fun16_5',rand(4,1),goal,abs(goal),a,b,[],[],zeros(4,1))

Fun16_5函数（可修改）

function F=Fun16_5(x)

F=[-100*x(1)-90*x(2)-80*x(3)-70*x(4);3*x(2)+2*x(4)];

end

Zg·ln

关注

5
点赞
踩
10

收藏

觉得还不错? 一键收藏
打赏
1
评论
如何构造求解规划模型——数学建模（MATLAB）

规划模型主要分为线性规划、整数规划、混合整数规划、非线性规划、多目标规划。接下来会一一讲述各类使用方法：可以直接用！！！
复制链接

扫一扫

专栏目录