机会约束程序的凸近似（学习用）（Matlab代码实现）

最新推荐文章于 2024-05-06 11:10:00 发布

然哥爱编程

最新推荐文章于 2024-05-06 11:10:00 发布

阅读量186

点赞数

文章标签：概率论算法机器学习

本文链接：https://blog.csdn.net/m0_73907476/article/details/133787729

版权

💥1 概述

📚2 运行结果

🎉3 参考文献

👨‍💻4 Matlab代码

💥1 概述

我们考虑一个机会约束问题，其中寻求在满足随机扰动凸约束系统的概率接近一的解上最小化凸目标。这个问题可能碰巧在计算上很难解决；我们的目标是建立其可计算处理的近似，即在机会约束问题中包含可行集的有效可解确定性优化程序。我们构造了相应的机会约束问题的一类一般的凸保守近似。

此外，在假设扰动中的约束是仿射的，并且扰动向量中的条目彼此独立于其他随机变量的情况下，本文通过建立机会约束问题的大偏差型近似，称为“Bernstein近似”。这种近似是凸的并且是有效可解的。本文提出了一种基于模拟的方案来界定机会约束问题中的最优值，并报告了旨在比较Bernstein和众所周知的情景近似方法的数值实验。最后，将构造扩展到模糊机会约束问题的情况，其中随机扰动与已知属于给定凸紧集而不是确切已知的分布的集合无关，而机会约束应该满足该集给定的每个分布。

📚2 运行结果

主函数部分代码：

%% Portfolio Problem

% Problem Value at Risk Minimation Problem
%  min J(x)=t
%  s.t. Prob{Sum(r_j x_j) < t} <= epsilon
%       Sum(x_j) = 1, x_j > 0 
clear all; close all;

d = 20;  % Number of assets

scenario  = 1;  % Conduct Scenario Appoach
bernstein = 0;  % Conduct Bernstein Approx

%% 1. Nonlinear Optimization with Scenario Approach
% Sample N samples and solve
% %  min J(x)=t
%  s.t. Prob{Sum(r_j x_j) < t} <= epsilon
%       Sum(x_j) = 1, x_j > 0 

if scenario
    eps_array  = 0.01:0.005:0.05;

    ncase      = length(eps_array);
    var_array  = zeros(1,ncase);
    vr_array   = zeros(1,ncase);

    rng(1);

    for ii = 1:ncase
        epsilon = eps_array(ii);
        beta    = 1 -epsilon;

        % Number of required samples
        eps_N = 0.01;
        N = ceil(2/(1-beta)*log(1/eps_N) + 2*d + 2*d/(1-beta)*log(2/(1-beta)));

        % create random variables
        % r = (ones(d,N) + 0.06*jvec) + 0.04*(jvec.*xi);  % d X N  non-affine
        r = weight(d,N);  % d X N affine

        % Test
        x0 = 1/d * ones(d,1);
        y0 = [x0; 2];

        % Linear Programming
        fobj = [zeros(1,d),-1];
        A   = [-r', ones(N,1)];  %N x (d+1) t - rx <= 0  
        b   = zeros(N,1);
        Aeq = [ones(1,d),0];
        beq = 1;
        lb  = zeros(d+1,1);
        ub  = [ones(d,1);2];

        % Solutions
        ysol = linprog(fobj,A,b,Aeq,beq,lb,ub);
        xsol = ysol(1:end-1);
        varsol = ysol(end);

        % calculate violation rate
        Nv = 10^6;
        vr = Vrate(xsol,varsol,Nv);

        % Print Results
    %     disp(['eps:',num2str(epsilon)]);
    %     disp(['    VaR:',num2str(tsol)]);
    %     disp(['    Violation Rate:',num2str(vr)]);

        var_array(ii) = varsol;
        vr_array(ii) = vr;
    end

    figure(1);
    hold on; grid on;
    title('Sampling Method');
    xlabel('eps')
    ylabel('VaR')
    plot(eps_array,var_array,'LineWidth',4);
    set(gca,'FontSize',14);

    figure(2);
    hold on; grid on;
    title('Sampling Method');
    xlabel('eps')
    ylabel('Violation Rate')
    plot(eps_array,vr_array,'LineWidth',4);
    set(gca,'FontSize',14);

end