Esta nota técnica mostra como utilizar a função solve.QP do package quadprog do R para solução de problemas de programação quadrática com restrições. Essa função implementa o método dual descrito em Goldfarb and Idnani (1982, 1983). Para detalhes específicos e outras opções, consulte a documentação do package.
\(\def\X{{\mathbf X}} \def\x{{\mathbf x}} \def\d{{\mathbf d}} \def\D{{\mathbf D}} \def\A{{\mathbf A}} \def\b{{\mathbf b}} \def\vecR{{\vec{R}}} \def\vectheta{{\vec{\theta}}} \def\vecmu{{\vec{\mu}}} \def\bfa{{\boldsymbol{a}}} \def\bfS{{\boldsymbol{\Sigma}}} \def\bfV{{\mathbf{V}}} \def\var{{\rm Var}} \def\cov{{\rm Cov}} \def\CV{{\rm cv}} \def\E{{\rm E}} \def\B{{\rm b}} \DeclareMathOperator*{\min}{min} \displaystyle \min_{\x}\;\; -\d^T\x+\frac{1}{2}\,\x^T\D \x\;\;\\ \text{sujeito a}\;\;\; \A^T\x \ge\;\mbox{e/ou}\; = \b\;\;\;\text{ (primeiras k restrições podem ser igualdades)}\)
na especificação do problema na função solve.QP há a possibilidade das primeiras k restrições serem igualdades em lugar de desigualdades do tipo \(\ge\). Basta especificar meq=k na função para indicar a existência de igualdades. Caso não seja especificado um valor para meq, será assumido meq=0.
uma restrição do tipo \(a_{11} x_1 + a_{12} x_2 \le b_1\) pode ser codificada como \(-a_{11} x_1 - a_{12} x_2 \ge -b_1\)
no problema as letras minúsculas \(\d\), \(\x\) e \(\b\) indicam vetores e as letras maiúculas \(\D\) e \(\A\) indicam matrizes.
na função, as variáveis associadas aos vetores e matrizes da formu lação geral, seguem uma regra, onde a terminação “vec” ou “mat” é incorporada a letra na fórmula (ex. Dmat, Amat, bvec, dvec)
para um problema de maximização (max) multiplique por -1 a função objetivo.
Suponha que deseja solucionar o problema envolvendo o modelo média-variância no contexto da otimização de carteiras de investimento:
\(\DeclareMathOperator*{\min}{min} \displaystyle \min_{\vectheta}\; \vectheta\,^{\small T}\,\bfV\,\vectheta\) sujeito a \(\vectheta\,^{\small T}\,\vecmu=\mu_0\;\;\;\sum_{i=1}^m \theta_i =1\;\;\;\mbox{e}\;\;\;\theta_i\ge 0\), \(i=1,\ldots,m\).
onde \(\vectheta=\left(\begin{array}{c} \theta_1\\ \theta_2\\ \theta_3\\ \end{array}\right)\) \(\vecmu=\left(\begin{array}{c} 0{,}12\\ 0{,}09\\ 0{,}08\\ \end{array} \right)\) \(\bfV=\left(\begin{array}{ccc} 4&-1{,}4&-1\\ -1{,}4&1&0{,}5\\ -1&0{,}5&1\\ \end{array}\right)\;\;\;\text{e}\;\;\;\mu_0=0{,}10\)
Na notação do problema, já na forma matricial, as restrições \(\vectheta\,^{\small T}\,\vecmu=\mu_0\;\;\;\sum_{i=1}^m \theta_i =1\;\;\;\mbox{e}\;\;\;\theta_i\ge 0\), \(i=1,\ldots,m\) poderiam ser representadas por:
\(\begin{array}{ccccccc} \theta_1 0{,}12&+&\theta_2 0,09&+&\theta_3 0,08&=&0{,}10\\ \theta_1&+&\theta_2&+&\theta_3&=&1\\ \theta_1&&&&&\ge&0\\ &&\theta_2&&&\ge&0\\ &&&&\theta_3&\ge&0\\ \end{array}\)
que na forma matricial seriam representadas por
Nessa última expressão, já usando a notação da função solve.QP, temos
\(\A^T = \left(\begin{array}{ccc} 0{,}12&0{,}09&0{,}08\\ 1&1&1\\ 1&0&0\\ 0&1&0\\ 0&0&1\\ \end{array}\right)\;\; \mbox{e}\;\;\) \(\b= \left( \begin{array}{c} 0{,}10\\ 1\\ 0\\ 0\\ 0\\ \end{array}\right)\;\;\;\mbox{onde}\;\;\;\) \(\A=\left(\begin{array}{ccccc} 0{,}12&1&1&0&0\\ 0{,}09&1&0&1&0\\ 0{,}08&1&0&0&1\\ \end{array}\right)\)
Note que no problema temos 2 restrições do tipo “=” e 3 restrições do tipo “\(\ge\)”.
Com a formatação das restrições no formato matricial realizada nos parágrafos anteriores, podemos reescrever o problema de minimização como
Agora basta converter esse modelo para a notação exigida pela função solve.QP faça as seguintes substituições:
\(\x=\vectheta\) e \(\D=2 \bfV\)
Defina agora o vetor
\(\left(\begin{array}{c} 0\\ 0\\ 0\\ \end{array}\right)\)
Com essas definições chegamos ao problema equivalente a ser resolvido em \(\x\)
\(\DeclareMathOperator*{\min}{min} \displaystyle \min_{\x}\; -\d^T \x \frac{1}{2}\x\,^{\small T}\,\D\,\x\\ \mbox{sujeito a}\\ \A^T\vectheta \left[ \begin{array}{c} =\\ =\\ \ge\\ \ge\\ \ge\\ \end{array}\right] \b\)
Passando agora à codificação:
## definição dos dados do problema
m<-c(0.12,0.09,0.08)
V<-matrix(c(4,-1.4,-1,-1.4,1,0.5,-1,0.5,1),nrow=3,ncol=3)
V## [,1] [,2] [,3]
## [1,] 4.0 -1.4 -1.0
## [2,] -1.4 1.0 0.5
## [3,] -1.0 0.5 1.0A<-matrix(c(0.12,0.09,0.08,1,1,1,1,0,0,0,1,0,0,0,1),nrow=3,ncol=5)
A## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.12 1 1 0 0
## [2,] 0.09 1 0 1 0
## [3,] 0.08 1 0 0 1## ou na forma equivalente
A<-matrix(c(m,1,1,1,1,0,0,0,1,0,0,0,1),nrow=3,ncol=5)
b<-c(0.10,1,0,0,0)
d<-c(0,0,0)
## Codificação usando a notação da função solve.QP
Dmat<-2*V
dvec<-d
Amat<-A
bvec<-b
numigual<-2 ## número de igualdades nas restrições
###Solução completa
require(quadprog) ## carregando o package (instale se necessário)## Warning: package 'quadprog' was built under R version 3.1.3solve.QP(Dmat, dvec, Amat, bvec, meq=numigual)## $solution
## [1] 0.35135135 0.59459459 0.05405405
##
## $value
## [1] 0.2594595
##
## $unconstrained.solution
## [1] 0 0 0
##
## $iterations
## [1] 3 0
##
## $Lagrangian
## [1] 25.945946 2.075676 0.000000 0.000000 0.000000
##
## $iact
## [1] 1 2## ou para recuperar somente o valor das variáveis na minimização
solve.QP(Dmat, dvec, Amat, bvec, meq=numigual)$solution## [1] 0.35135135 0.59459459 0.05405405## o valor mínimo obtido é recuperado por
solve.QP(Dmat, dvec, Amat, bvec, meq=numigual)$value## [1] 0.2594595Assim podemos concluir que o valor mínimo da função objetivo (atendendo as restrições) é 0,2225, obtido com