Programmation quadratique en python à partir des données d'octave

Question

Je suis en train de convertir des parties d'un programme MATLAB en Python et Octave.

J'utilise Octave pour générer deux matrices, puis j'importe ces matrices dans python en utilisant oct2py . La racine de mon problème est constituée par ces lignes dans MATLAB ( H_combined y f_combined ci-dessous)

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

Actuellement, je cherche une solution en Python ou Octave qui produirait un résultat similaire, mais en vain.

J'ai tenté d'utiliser quadprog de Octave optim Cependant, j'obtiens le résultat suivant 120, 1, 1, 1, ..., 1 sur x_ridge_combined au lieu d'un assortiment de valeurs flottantes, comme je m'y attendais. J'ai vérifié que H_combined y f_combined sont exactement les mêmes que lorsqu'ils sont exécutés dans MATLAB, mais je suppose qu'il est possible d'obtenir des informations plus détaillées. quadprog dans Octave ne fonctionne pas de la même manière.

Après avoir essayé une approche Octave, je tente d'importer les valeurs dans Python pour essayer d'utiliser la fonction quadprog paquete.

Essai quadprog ,

print(quadprog.solve_qp(H,f))

donne l'erreur

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

Les types et formes de H y f sont les suivants :

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

Quelqu'un sait-il pourquoi j'obtiens ces erreurs ? Ou d'autres suggestions sur la façon de traduire cette ligne à partir de MATLAB ?

Answer 1

Essayez d'utiliser cvxopt_quadprog . L'auteur prétend qu'il imite MATLAB quadprog et il devrait accepter les arguments de la même manière :

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
    """

L'erreur est probablement due au fait que votre f est une matrice [1x123], alors qu'il devrait s'agir d'un vecteur de longueur [123]. Vous pouvez essayer de la remodeler :

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

Answer 2

Oui, mais le problème de cvxopt_quadprog est qu'il est considérablement plus lent pour les grandes optimisations itératives de séries temporelles car il vérifie à chaque fois si le problème est PSD. C'est pourquoi j'espérais utiliser quad_prog, qui s'est avéré beaucoup plus rapide. Ref : https://github.com/stephane-caron/qpsolvers

Answer 3

Bien qu'il s'agisse d'un projet un peu en dehors du champ d'application, je souhaite apporter au projet NLopt en jeu. Comme son acronyme l'indique, il s'attaque à l'optimisation non linéaire, mais avec de nombreuses possibilités. algorithmes globaux/locaux, sans dérivées/avec dérivées explicites . La raison pour laquelle je souhaite le mentionner est qu'il dispose d'une interface pour MATLAB, Octave + Python (et C/C++,...). Il est donc très facile de reproduire les solutions dans différents langages (c'est la raison pour laquelle je l'ai découvert) ; de plus, les algorithmes sont en fait plus rapides que les algorithmes natifs de MATLAB (c'est ma propre expérience).

Pour votre problème, je choisirais BOBYQA (optimisation limitée par optimisation quadratique ou SLSQP (programmation quadratique séquentielle par moindres carrés). Cependant, vous devrez écrire une fonction de coût plutôt que de remettre des matrices

L'installation est facile via pip

pip install nlopt

faire une petite vérification

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

un code rapide sur la façon d'utiliser l'optimisation :

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)

Dans MATLAB, il suffit d'ajouter les DLL au chemin d'accès. J'ai écrit un petit wrapper pour BOBYQA afin d'imiter l'interface de MATLAB (au cas où vous voudriez le tester dans les deux langues =P -- faites-le moi savoir, je l'utilise plus souvent dans MATLAB... comme le wrapper le montre probablement^^) :

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