用Octave的数据在Python中进行二次规划

Question

我正在将MATLAB程序的部分转换为Python和Octave。

我正在使用Octave生成两个矩阵，然后使用oct2py将这些矩阵导入Python。 我问题的根源就是下面这些MATLAB代码中的这两行（H_combined和f_combined）。

handles.options =optimset('algorithm','interior-point-convex','Display','off','TolFun',1e-15,'TolX',1e-10,'MaxFunEvals', 1E5);

handles.x_ridge_combined = quadprog(H_combined, f_combined, [], [], [], [], handles.lb_re, handles.ub_re, handles.x_re_0, handles.options);

目前，我正在寻找Python或Octave的解决方案，以产生类似的输出，但未能找到合适的。

我尝试使用Octave的optim中的quadprog，但是在x_ridge_combined上得到了一个输出为120, 1, 1, 1, ..., 1的结果，而不是我所期望的一系列浮点值。 我已经验证了H_combined和f_combined与在MATLAB中运行时完全相同，但我想Octave中的quadprog并不起作用。

在尝试了Octave的方法之后，我尝试将这些值导入Python，并尝试使用quadprog包。

尝试使用quadprog，

print(quadprog.solve_qp(H,f))

产生了错误

ValueError: Buffer has wrong number of dimensions (expected 1, got 2)

H和f的类型和形状如下：

<class 'numpy.ndarray'> #H
(123, 123)
<class 'numpy.ndarray'> #f
(1, 123)

谁知道为什么我会收到这些错误？或者对如何从MATLAB中翻译那一行有其他的建议吗？

- nrobins1

---

快速阅读后，我认为您最好提供更多信息。您的matlab原始代码与您合并的Python/Octave替代品之间的分歧点似乎是quadprog。即使使用完全相同的输入，Python/Octave的quadprog也无法产生与matlab相同的输出。这是正确的吗？ 您能否发布您正在使用的输入？ 如果这些是大型数据集，请您生成一个可以产生相同错误的较小示例吗？ 我知道这需要一些工作... - sancho.s ReinstateMonicaCellio

Matlab中的quadprog命令和Octave中的完全一样吗？ - sancho.s ReinstateMonicaCellio

尝试使用Matlab和Octave中的文档示例来检查结果。 https://octave.sourceforge.io/optim/function/quadprog.html - sancho.s ReinstateMonicaCellio

请查看这个关于Matlab和Octave的区别的帖子，虽然我不确定它是否适用于您的情况 https://lists.gnu.org/archive/html/help-octave/2015-09/msg00120.html - sancho.s ReinstateMonicaCellio

或者检查这些其他示例以进行比较https://www.mathworks.com/help/optim/ug/quadprog.html。 - sancho.s ReinstateMonicaCellio

Answer 1

尝试使用cvxopt_quadprog。作者声称它模仿了MATLAB quadprog，并且应该接受相同的参数：

def quadprog(H, f, L=None, k=None, Aeq=None, beq=None, lb=None, ub=None):
    """
    Input: Numpy arrays, the format follows MATLAB quadprog function: https://www.mathworks.com/help/optim/ug/quadprog.html
    Output: Numpy array of the solution
    """

很可能出错是因为您的f矩阵是[1x123]，而应该是长度为[123]的向量。您可以尝试重新整形它：

f = f.reshape(f.shape[1])

- igrinis

---

我同意“reshape”解决方案。应该适用于原始代码 =) - max

Answer 2

是的，虽然cvxopt_quadprog存在一个问题，就是对于大型迭代时间序列优化，它会检查每次问题是否为PSD，所以速度相对较慢。因此，我希望能够利用已被证明更快的quad_prog。请参见: https://github.com/stephane-caron/qpsolvers

- Panda_User

Answer 3

虽然有点偏题，但我想提一下NLopt项目。正如首字母缩写所示，它处理非线性优化问题，但具有丰富的全局/局部、无导数/显式导数算法。我想提到它的原因是，它有MATLAB、Octave+Python（以及C/C++）的接口。因此，它很容易在不同的语言中重现解决方案（这就是我发现它的原因）；此外，算法实际上比MATLAB本地算法更快（这是我的个人经验）。

对于您的问题，我建议使用BOBYQA（二次规划有界优化）或SLSQP（顺序最小二乘二次规划）。但是，您需要编写成本函数而不是交出矩阵。

通过pip安装非常容易

pip install nlopt

做一个小检查

import nlopt
# run quick test. Look for "Passed: optimizer interface test"
nlopt.test.test_nlopt()

关于如何使用优化的一些快速代码：

import numpy as np
import nlopt

obj = nlopt.opt(nlopt.LN_BOBYQA,5)
obj.set_min_objective(fnc)

obj.set_lower_bounds(lb)
obj.set_upper_bounds(ub)

def fnc(x, grad):

        """
        The return value should be the value of the function at the point x, 
        where x is a NumPy array of length n of the optimization parameters 
        (the same as the dimension passed to the constructor).

        In addition, if the argument grad is not empty, i.e. grad.size>0, then 
        grad is a NumPy array of length n which should (upon return) be set to 
        the gradient of the function with respect to the optimization parameters 
        at x. That is, grad[i] should upon return contain the partial derivative , 
        for , if grad is non-empty.
        """
        H = np.eye(len(x)) # extampe matrix

        cost = 0.5*x.transpose().dot( H.dot(x) )
        return float(cost) # make sure it is a number

xopt = obj.optimize(x0)

在MATLAB中，您只需要将DLL添加到路径中即可。我编写了一个简短的包装器来模仿MATLAB的界面（如果您想在两种语言中进行比较，请告诉我，我更经常在MATLAB中使用它...因为包装器可能会显示^^）：

function [x_opt, fval, exitflag] = BOBYQA(fnc,x0,lb,ub, varargin)
% performes a constrained, derivative-free local optimization
%
% --- Syntax:
% x_opt = BOBYQA(fnc,x0,lb,ub)
% x_opt = BOBYQA(...,'MaxEval',10)
% x_opt = BOBYQA(...,'MaxTime',5)
% [x_opt, fval] = BOBYQA(...)
% [x_opt, fval, exitflag] = BOBYQA(...)
% 
% --- Description:
% x_opt = BOBYQA(fnc,x0,lb,ub)  takes a function handle 'func', an initial 
%               value 'x0' and lower and upper boundary constraints 'lb' 
%               and 'ub' as input. Performes a constrained local 
%               optimization using the algorithm BOBYQA from Powell 
%               http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2009_06.pdf.
%               Returns the optimal value 'x_opt'.
% x_opt = BOBYQA(...,'MaxEval',10)optional input parameter that defines the
%               maximum number of evaluations.
% x_opt = BOBYQA(...,'MaxTime',5) optional input parameter that defines the
%               maximum allowed time in seconds for the optimization. This 
%               is a soft constraint and may be (slightly) broken.
% [x_opt, fval] = BOBYQA(...) seconds return value is the optimal function 
%               value.
% [x_opt, fval, exitflag] = BOBYQA(...) third return value is the exitflag,
%               see function NLoptExitFlag().
% 
% ------------------------------------------------------------------- 2017

% NLOPT_LN_BOBYQA

% --- parse input
IN = inputParser;
addParameter(IN,'MaxEval',10000, @(x)validateattributes(x,{'numeric'},{'positive'}));
addParameter(IN,'MaxTime',60, @(x)validateattributes(x,{'numeric'},{'positive'}));
parse(IN,varargin{:});

% generic success code: +1
%      stopval reached: +2
%         ftol reached: +3
%         xtol reached: +4
%      maxeval reached: +5
%      maxtime reached: +6
% generic failure code: -1
%    invalid arguments: -2
%        out of memory: -3
%     roundoff-limited: -4

    % set options
    opt = struct();
    opt.min_objective = fnc;
    opt.lower_bounds = lb;
    opt.upper_bounds = ub;



    % stopping criteria
    opt.maxtime = IN.Results.MaxTime; % s  % status = +6
%     opt.fc_tol = FncOpt.STOP_FNC_TOL*ones(size(ParInit)); % +3
%     opt.xtol_rel = FncOpt.STOP_XTOL_REL; % +4
%     opt.xtol_abs = FncOpt.STOP_XTOL_ABS*ones(size(ParInit)); % +4
    opt.maxeval = IN.Results.MaxEval; % status = +5

    % call function
    opt.algorithm = 34;% eval('NLOPT_LN_BOBYQA');

    t_start = tic;
    [x_opt, fval, exitflag] = nlopt_optimize(opt,x0);
    dt = toc(t_start);
    fprintf('BOBYQA took %.5f seconds | exitflag: %d (%s)\n',dt,exitflag,NLoptExitFlag(exitflag))
end

function txt = NLoptExitFlag(exitflag)
% generic success code: +1
%      stopval reached: +2
%         ftol reached: +3
%         xtol reached: +4
%      maxeval reached: +5
%      maxtime reached: +6
% generic failure code: -1
%    invalid arguments: -2
%        out of memory: -3
%     roundoff-limited: -4

switch exitflag
    case 1
        txt = 'generic success code';
    case 2
        txt = 'stopval reached';
    case 3
        txt = 'ftol reached';
    case 4
        txt = 'xtol reached';
    case 5
        txt = 'maxeval reached';
    case 6
        txt = 'maxtime reached';
    case -1
        txt = 'generic failure code';
    case -2
        txt = 'invalid arguments';
    case -3
        txt = 'out of memory';
    case -4
        txt = 'roundoff-limited';
    otherwise
        txt = 'undefined exitflag!';
end
end

- max