数学建模中的规划模型总结：MATLAB求解方法

function lp_
% MATLAB求解线性规划问题，linprog
%
c = -[4000;3000]; % 目标函数系数
A = [2 1;1 1;0 1];
b = [10;8;7];
Aeq = [];
beq = [];
lb = [0 0]; % 变量下界
x = linprog(c,A,b,Aeq,beq,lb)
fval = -c'*x

% 返回的fval是负的，
% 即使问题是正的，
% 在matlab内部，prob2struct将最大化问题转化为目标函数负的最小化问题

x1 = optimvar('x1','LowerBound',0,'UpperBound',inf);
x2 = optimvar('x2','LowerBound',0,'UpperBound',7);
prob = optimproblem('Objective',4000*x1 + 3000*x2,'ObjectiveSense','max');

prob.Constraints.c1 = 2*x1 + x2 <= 10;
prob.Constraints.c2 = x1 + x2 <= 8;

problem = prob2struct(prob);
[sol,fval,~,output] = linprog(problem)

不包含非线性约束
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB)
包含非线性约束
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)
自己定义目标函数和非线性约束
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(@(X)MYOBJ(X),X0,A,B,Aeq,Beq,LB,UB,@(X)MYCON(X))性约束函数
% 目标函数
function F = MYOBJ(X)
F = ......
% 非线性约束函数
function [G,Heq] = MYCON(X)
G = ..... % 非线性不等式约束条件,G<=0
Heq = ..... % 非线性等式约束条件,Heq == 0

主函数
clear
n = 5;
x1 = [n;nan*ones(n,1)];
x2 = [0:n]';
x = [x1 x2 x2 x2]

[p f] = dynprog(x,'DecisF1','SubObjF1','TransF1','ObjF1')

disp('p的第三列是决策变量值')
决策函数

function u = DecisF1(k,x)
% 决策函数
if k ==4
    u =x;
else
    u = 0:x;
end
end
状态转移方程

function y = TransF1(k,x,u)
% 状态转移方程
y = x-u;
end
阶段k的指标函数
function v = SubObjF1(k,x,u)
% 阶段k的指标函数
a = [0.2 0.3 0.5 0.9];
A = [8 7 6 3];
v = A(k)*(1-exp(-a(k)*u));
v = -v;
end
基本方程中的函数
function y = ObjF1(v,f)
% 基本方程中的函数
y = v+f;
% y = -y;
end
动态规划dynprog函数
function [p_opt, fval] = dynprog(x,DecisFun, SubObjFun, TransFun, ObjFun)

% 动态规划逆序算法

% x是状态变量，一列代表一个阶段状态；
% DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量；
% SubObjFun(k,x,u)是阶段指标函数;
% TransFun(k,x,u)是状态转移函数，其中x是阶段k的某状态变量，u是相应的决策变量；
% ObjFun(v,f)是第k阶段至最后阶段指标函数，当ObjFun(v,f) = v+f时，输入的ObjFun可以省略；
% fval是一个列向量，各元素分别表示p_opt各最优策略组对应始端状态x的最优函数值

function path=AStar(obstacle,map)
path=[];        %用于存储路径
open=[];        %OpenList,考虑的节点
close=[];       %CloseList，不考虑的节点
findFlag=false; %findFlag用于判断while循环是否结束

%% 1.将起始点放在Openlist中

% open：[节点坐标，代价值F=G+H,代价值G,父节点坐标]
open =[map.start(1), map.start(2) ,...
    0+h(map.start,map.goal) ,...
    0 , ...
    map.start(1) , map.start(2)];

%与当前节点的坐标差值（前两列）
%与当前节点的距离值（最后一列）

NEXT = [-1,1,14;...
    0,1,10;...
    1,1,14;...
    -1,0,10;...
    1,0,10;...
    -1,-1,14;...
    0,-1,10;...
    1,-1,14];

%% 2. 循环以下过程

while ~findFlag

    %--------------------首先判断是否达到目标点，或无路径-----
    if isempty(open(:,1))
        disp('没有目标路径');
        return;
    end

    %判断目标点是否出现在open列表中
    [isopenFlag,Id]=isopen(map.goal,open);
    if isopenFlag
        disp('找到目标');
        close = [open(Id,:);close];
        findFlag=true;
        break;
    end

    %按照Openlist中的第三列（代价函数F）进行排序，查找F值最小的节点
    [Y,I] = sort(open(:,3)); %对OpenList中第三列排序
    open=open(I,:);%open中第一行节点是F值最小的

    %将F值最小的节点(即open中第一行节点)，放到close第一行(close是不断积压的)，作为当前节点
    close = [open(1,:);close];
    current = open(1,:);
    open(1,:)=[];%因为已经从open中移除了，所以第一列需要为空

    %对当前节点周围的8个相邻节点进行计算，算法的主体：------------------------
    for in=1:length(NEXT(:,1))

        m=[current(1,1)+NEXT(in,1) , current(1,2)+NEXT(in,2) , 0 , 0 , 0 ,0];
        m(4)=current(1,4)+NEXT(in,3); % m(4)  相邻节点G值
        m(3)=m(4)+h(m(1:2),map.goal);% m(3)  相邻节点F值

        if isObstacle(m,obstacle)
            continue;
        end

        [flag,targetInd]=findIndex(m,open,close);

        if flag==1
            continue;
        elseif flag==2
            m(5:6)=[current(1,1),current(1,2)];
            open = [open;m];
        else
            if m(3) < open(targetInd,3)
                m(5:6)=[current(1,1),current(1,2)];
                open(targetInd,:) = m;
            end
        end
    end

    %绘制
    pause(0.01);
    fillPlot(close,[204,159,42]/255);
    hold on;
    fillPlot(open,[51,102,153]/255);
end

%% 3. 追溯路径
path=GetPath(close,map.start);
end

% 载入地图
% --------------------------------------------------------—————————
% PIC = imread(uigetfile('*.png'));         % 载入地图

PIC = imread('map2.png');                    % 载入地图
imshow(PIC);                                % 显示地图
hold on;

% 关键参数
% --------------------------------------------------------—————————
feasibleR = 51;                             % 可行区域的红色值【0~255】
NUM = 200;                                  % 离散点规模

% 一些样式设置，可以注释掉
% --------------------------------------------------------—————————
ax = gca;
ax.Parent.Units = 'normalized';
ax.Parent.Color = 'k';
ax.Parent.Position = [0.1 0.1 0.8 0.7];
ax.Units = 'normalized';
ax.Position = [0.1 0.001 0.8 0.998];
ax.Color = [51 51 51]/255;

% 选择起/始点
% --------------------------------------------------------—————————
StartPoint = ginput(1);
EndPoint = ginput(1);

scatter([StartPoint(1) EndPoint(1)],[StartPoint(2) EndPoint(2)],...
    'filled','SizeData',150,'MarkerFaceColor',[199 237 53]/255,'MarkerEdgeColor','w');
% 绘制出起始点
pause(1)

% 生成并绘制可行点
%__________________________________________________________________________
Height = size(PIC,1);
Width = size(PIC,2);
randPoints = floor([rand(NUM,1)*Height,...
    rand(NUM,1)*Width])+1;                  % 空间随机撒点
allPoints = [StartPoint;EndPoint];          % 初始化可行点列表
for i=1:NUM
    if(PIC(randPoints(i,2),randPoints(i,1),1) == feasibleR  )
        allPoints = [allPoints;double(randPoints(i,:))];
    end
end% 得到所有可行点列表
for i = 1:size(allPoints,1)
    scatter(allPoints(i,1),allPoints(i,2),'filled','SizeData',20,'MarkerFaceColor',[199 237 53]/255);
    pause(.0005)
end

% 初始化无向图邻接矩阵
%__________________________________________________________________________
AdjacencyMatrix = genAM(PIC,allPoints);

% 对矩阵AdjacencyMatrix使用【Dijkstra最短路径算法】；
%__________________________________________________________________________