Unity开发者联盟收藏本站
发帖教程新人必读切换到宽版

 找回密码
 立即注册
Unity开发者联盟»论坛 › 3D Answers 技术讨论 › 技术讨论 Technology Answers › 凸优化课程设计1
返回列表 发新帖
查看: 412|回复: 0

凸优化课程设计1

[复制链接]
发表于 2022-6-28 18:37 | 显示全部楼层 |阅读模式
Project Tasks:
1、运用编程语言(推荐使用Matlab)，设计完成一个通用的单纯形算法，完成给定线性规划问题的求解，并与Matlab自带的函数linprog对比算法结果和运行效率。
2、运用编程语言(推荐使用Matlab)，设计完成一个通用的分枝定界算法，完成给定整数规划问题的求解，并与Matlab自带的函数 intlinprog对比算法结果和运行效率。
Requirements:
1、用单纯形法求解下述线性规划问题(要求写出步骤)，并与编程实现的通用单纯形算法对比，验证所完成的单纯形算法的正确性


2、用分枝定界算法求解下述整数规划问题(要求写出步骤)，并与编程实现的通用分枝定界算法算法对比，验证所完成的分枝定界算法的正确性


3、用验证后的分枝定界算法，完成下述5个城市的TSP问题，其中5个城市间距离如下表所示


1. 单纯形算法

1. 问题建模



2. 单纯形算法流程



3. 辅助线性问题流程



4. 算法对比

4.1 效率对比
       运行linprog函数用时如图1-1所示：0.148s
       运行Laux函数如图1-2所示：0.033s
       分析：simplex函数在Laux函数里使用。Laux负责判断原问题是否有最优解，有则给出一个可行解送给simplex函数。simplex函数未完善，不能解决有等式约束的问题并且未解决死循环自动退出的问题。



图1-1



图1-2

4.2 结果对比
       图1-3和图1-4分别为运行linprog函数与Laux函数的过程。图1-5和图1-6分别为linprog函数与Laux函数的输出结果，可见结果一致。



图1-3



图1-4



图1-5



图1-6

2. 分支定界算法

1. 问题建模



2. 分支定界算法流程



迭代方法：
1. 分，在线性规划问题A的最优解中任选个不符合整数条件的变量xj，其值为bj，构造两个约束条件，xj=floor(bj)、xj=ceil(bj)，分别加到问题A中，形成两个问题。不考虑整数条件求解这两个问题。即分
2. 定界。对每个后继问题表明其求解的结果，与其他问题进较，将最优标函数值最者（不包括问题A）作为新的下界，在已符合整数条件的各分中，找出标函数值最者作为新的上界。
3. 剪枝，将标函数值不在上界、下界中的分剪去
4. 重复1 2 3，直到得到最优解
注意事项：在对两分进分解时，优先选择标函数值最的分进分解。
3. 算法对比

3.1 效率对比
       运行intlinprog函数用时如图2-1所示：0.152s
       运行LP函数如图2-2所示：0.293s
       分析：因为单纯形算法写得不完善（未解决有等式约束时的问题），所以使用linprog函数求最优值。因为分定界法的求解过程与叉树结构相同，所以利叉树这种数据结构来存储每个线性规划问题的解及对应的标函数值。



图2-1



图2-2

3.1 结果对比
       图2-3和图2-4分别为运行intlinprog函数与LP函数的过程。图2-5和图2-6分别为intlinprog函数与LP函数的输出结果，可见结果一致。



图2-3



图2-4



图2-5



图2-6

3. TSP问题

1. 问题建模



2. 程序代入

       分析：图3-1为使用LP函数解决TSP问题的过程，初定义的a矩阵对应 （意为每个起点只能有一个终点），末定义的矩阵coff对应 （意为每个终点只能有一个起点）。图3-2为输出结果，即取f的2、5、12、14、20列相加，这些列数分别对应x13、x21、x35、x42、x54。



图3-1



图3-2

附录

1.单纯形算法代码

Laux为主函数，调用simplex.m。
function [x,z,ST,res_case] = Laux(c,A,b)
% ------------------------------------------------------------------------------------------------------
% 先把线性规划转换为标准模型：1.目标函数求最小值；2.约束条件全部为小于等于不等式；3.决策变量非负.
% 未解决问题：1.未去冗余 2.未考虑等式约束
% 输入参数: f为目标函数系数, A为约束不等式方程组系数矩阵, b为约束不等式方程组常数项(标准线性规划)
% 输出参数: x最优解, z最优目标函数值, ST存储单纯形表数据, res_case=1表示有最优解, res_case=0表示有无界解
% ------------------------------------------------------------------------------------------------------
%% --------------去冗余---------------



%% ----------检测bi是否存在负值，有则改为正值--------(退化问题)
[m,n] = size(A);                   %m约束条件个数, n决策变量数(m≤n)
[m2,~] = size(b);                  %m2约束条件个数
ind_B = zeros(1,m);                 %初始化基变量索引
if any(b < 0)                      %进入退化问题
    coff_x0 = -1*ones(m,1);        %引入Laux
    new_A = [A,coff_x0];
    unit = eye(m);                 %松弛操作
    new_A = [new_A,unit];
    for i = 1:m                    %初始基可行解索引
        ind_B(i) = n+1+i;
    end
    [~,y] = min(b);                %y为负值最大的索引
    ind_N = setdiff(1:m2, y);      %非基变量的索引
    for i = 1:length(ind_N)        %把bi修改为正值
        new_A(ind_N(i),:) = new_A(ind_N(i),:) - new_A(y,:);
        b(ind_N(i)) = b(ind_N(i)) - b(y);
    end
    new_A(y,:) = -new_A(y,:);
    b(y) = -b(y);
    ind_B(1) = n+1;
    %%%%%%%%%%%%%%%找出原问题的初始基可行解%%%%%%%%%%%%%
    for i = 1:n
        if abs(new_A(y,n+1)-1.0)<0.001 && all(abs(new_A(ind_N,n+1)-0.0))<0.001  %判断x0是否为基变量,是就pivot，不是就跳出循环得到基变量
            theta = b./(new_A(:,i));
            [~,pivot_index] = min(theta(find(theta>0)));
            nopiv_index = setdiff(1:m, pivot_index);
            ind_B(pivot_index) = i;                                             %更新基可行解索引
            for j = 1:m-1                                                       %pivot操作开始
                middle_coffe1 = new_A(nopiv_index(j),i);
                middle_coffe2 = new_A(pivot_index,i);
                new_A(nopiv_index(j),:) = new_A(nopiv_index(j),:) - middle_coffe1/middle_coffe2* new_A(pivot_index,:);
                b(nopiv_index(j),:) = b(nopiv_index(j),:) - middle_coffe1/middle_coffe2* b(pivot_index,:);
            end
            new_A(pivot_index,:) = new_A(pivot_index,:)/middle_coffe2;
            b(pivot_index,:) = b(pivot_index,:)/middle_coffe2;                  %pivot操作结束
        else
            break
        end
    end
    res_case = 1;
    if abs(new_A(y,n+1)-1.0)<0.000001 && all(abs(new_A(ind_N,n+1)-0.0))<0.001
        res_case = 0;                                                           %pivot结束后，x0仍为基变量，原问题无最优解
    end
    newA = [new_A(:,1:n),new_A(:,n+2:n+1+m)];                                   %A除去x0
    for i = 1:m
        if ind_B(i)>n                                                           %索引ind_B出去x0
            ind_B(i) = ind_B(i)-1;
        end
    end
else                         %不进入退化问题
    unit = eye(m);
    newA= [A,unit];          %松弛操作
    for i = 1:m
        ind_B(i) = n+i;
    end
    res_case = 1;
end

%% ---------调出Laux的基可行解送入单纯形法-------
c = [c,zeros(1,m)];
if res_case == 1;
    [x,z,ST,res_case] = simplex(c,newA,b,ind_B);   
    x([n+1:n+m]) = [];
end      

function [x,z,ST,res_case] = simplex(c,A,b,ind_B)
% -------------------------------------------------------------------------------------------------
% 先把线性规划转换为标准模型：1.目标函数求最小值；2.约束条件全部为小于等于不等式；3.决策变量非负.
% 输入参数: c为目标函数系数, A为约束方程组系数矩阵, b为约束方程组常数项, ind_B为基变量索引
% 输出参数: x最优解, z最优目标函数值, ST存储单纯形表数据, res_case=1表示有最优解, res_case=0表示有无界解
% 未解决问题：1.未解决退化问题 2.若f存在常数项，需要计算后手动添加
% -------------------------------------------------------------------------------------------------
%% --------------main--------------
[m,n] = size(A);              %m约束条件个数, n决策变量数
ind_N = setdiff(1:n, ind_B);  %非基变量的索引
ST = [];
% 循环求解
while true
    x0 = zeros(n,1);
    x0(ind_B) = b;               %初始基可行解（全部b值）
    cB = c(ind_B);               %计算cB
    Sigma = zeros(1,n);
    Sigma(ind_N) = c(ind_N) - cB*A(:,ind_N);   %计算检验数
    [~, k] = min(Sigma);         %选出最小检验数, 确定进基变量索引k
    Theta = b ./ A(:,k);         %计算θ
    Theta(Theta<=0) = 10000;
    [~, q] = min(Theta);         %选出最小θ
    el = ind_B(q);               %确定出基变量索引el, 主元为A(q,k)
    vals = [cB',ind_B',b,A,Theta];
    vals = [vals; NaN, NaN, NaN, Sigma, NaN];
    ST = [ST; vals];
    if ~any(Sigma < 0)           %此基可行解为最优解, any存在某个>0        
        x = x0;
        z = c * x;
        res_case = 1;
        return
    end
    if all(A(:,k) <= 0)          %有无界解
        x = [];
        res_case = 0;
        break
    end
    % 换基
    ind_B(ind_B == el) = k;      %新的基变量索引
    ind_N = setdiff(1:n, ind_B); %非基变量索引
    % 更新A和b
    A(:,ind_N) = A(:,ind_B) \ A(:,ind_N);
    b = A(:,ind_B) \ b;
    A(:,ind_B) = eye(m,m);
end2.分支定界代码

LP为主函数，调用createBinTreeNode.m、createBinTreeNode.m、IsInRange.m
function [x,z] = LP(f, A, b, Aeq, beq, lbnd, ubnd)
% --------------------------------------------------------------------------------------------------
% 先把线性规划转换为标准模型：1.目标函数求最小值；2.约束条件中的不等式全部为小于等于不等式
% 未解决问题：
% 输入参数: f为目标函数系数, A为约束不等式方程组系数矩阵, b为约束不等式方程组常数项(标准线性规划),
% Aeq为等式方程组系数矩阵, beq为等式方程组常数项, lbnd为决策变量下限, ubnd为决策变量上限
% 输出参数: x最优解, z最优目标函数值
% --------------------------------------------------------------------------------------------------
% 求解模型：
% min z = f * x
% A * x <= b
% x >= 0,且为整数
% ---------------------------------------------------
clear global;

global result; % 存储所有整数解
global lowerBound; % 下界
global upperBound; % 上界
global count; % 用以判断是否为第一次分支

count = 1;

BinTree = createBinTreeNode({f, A, b, Aeq, beq, lbnd, ubnd});

if ~isempty(result)
    [fval,flag]=min(result(:,end)); % result中每一行对应一个整数解及对应的函数值
    Result=result(flag,:);
%     disp('该整数规划问题的最优解为：');
    x = Result(1,1:end-1);
%     disp(x);
%     disp('该整数规划问题的最优值为：');
    z = Result(1,end);
%     disp(z);
else
    disp('该整数规划问题无可行解');
end
% createBinTreeNode.m
% 构建二叉树，每一结点对应一解
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 = [];
        % 寻找非整数解分量
        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));
            % 进行比较，优先选择目标函数较大的进行分支
            [~, 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               
                % 创建左子树
                BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});
                % 创建右子树
                BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
            else
                % 创建右子树
                BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
                % 创建左子树
                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

% IsInRange.m
% 判断分支问题的解是否在上下界的范围中，若不在，剪去
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
% 打印输出上下界
function y = OutputLowerAndUpperBounds()

global lowerBound;
global upperBound;

disp("此时下界、上界分别为");
disp(['下界：',num2str(lowerBound)]);
if isempty(upperBound)
    disp('无上界')
else
    disp(['上界：',num2str(upperBound)]);
end
3.TSP代码

构建约束条件的矩阵代入分支定界函数中
f = [2,1,3,4,1,4,4,2,5,4,2,2,5,2,2,3,4,2,4,2];
a = [ones(1,4),zeros(1,16);zeros(1,4),ones(1,4),zeros(1,12);zeros(1,8),ones(1,4),zeros(1,8);
zeros(1,12),ones(1,4),zeros(1,4);zeros(1,16),ones(1,4)];
coff = zeros(5,20);
coff(1,[5,9,13,17])=1;
coff(2,[1,10,14,18])=1;
coff(3,[2,6,15,19])=1;
coff(4,[3,7,11,20])=1;
coff(5,[4,8,12,16])=1;
a = [a;coff];
b = ones(10,1);
lbnd = zeros(20,1);
[x,z] = LP(f, [], [], a, b, lbnd, [])

本帖子中包含更多资源

您需要 登录 才可以下载或查看，没有账号？立即注册

×
回复

使用道具 举报

返回列表 发新帖
懒得打字嘛，点击右侧快捷回复 【右侧内容，后台自定义】
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Unity开发者联盟 ( 粤ICP备20003399号 )

GMT+8, 2024-11-15 20:33 , Processed in 0.092279 second(s), 26 queries .

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表