#quadratic programming
#install.packages("quadprog")
library(quadprog)

#min a - 2b + 4c + a^2 + 2b^2 + 3c^2 + ac
#subject to:
#3a + 4b - 2c <= 10
#-3a + 2b + c >= 2
#2a + 3b + 4c = 5
#0 <= a <= 5
#1 <= b <= 5
#0 <= c <= 5

#Package quadprog optimizes the following function
#min(-d^T b + 1/2 b^T D b)
#subject to:
#A^T b >= b_0

#rewrite optimization
#min a - 2b + 4c + a^2 + 2b^2 + 3c^2 + ac
#min -(-a+2b-4c) + 1/2(2a^2 + 4b^2 + 6c^2 +2ac)

#The first few constraints can be equalities (=), and the following constraints are all >=.

#Rewriting and reordering our constraints gives us
#2a + 3b + 4c = 5
#-3a - 4b + 2c >= -10
#-3a + 2b + c >= 2
#a >= 0
#-a >= -5
#b >= 1
#-b >= -5
#c >= 0
#-c >= -5

#Also rewrite the optimization function:

#The function call is defined as
#solve.QP(Dmat, dvec, Amat, bvec, meq)
#dvec = d
#Dmat = D
#bvec = b_0
#Amat = A
#meq is the number of equalities in constraints

#optimization variables
#min -(-a+2b-4c) + 1/2(2a^2 + 4b^2 + 6c^2 +2ac)
#we split 2ac to ac + ca!
dvec <- c(-1, 2, -4)
Dmat <- matrix(c(2, 0, 1, 
                 0, 4, 0,
                 1, 0, 6), ncol = 3, byrow = TRUE)

#constraints
Amat <- matrix(c(2, 3, 4,
                 -3, -4, 2,
                 -3, 2, 1,
                 1, 0, 0,
                 -1, 0, 0,
                 0, 1, 0,
                 0, -1, 0,
                 0, 0, 1,
                 0, 0, -1
), ncol = 3, byrow = TRUE)
Amat <- t(Amat)
bvec <- c(5, -10, 2, 0, -5, 1, -5, 0, -5)


dvec
Dmat

Amat
bvec

res <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
res

#manual check
a <- res$solution[1]
b <- res$solution[2]
c <- res$solution[3]
a - 2*b + 4*c +a*a + 2*b*b+3*c*c+a*c

#analytical solution
#a <- 499/1718
#b <- 1214/859
#c <- 77/1718
#a - 2*b + 4*c +a*a + 2*b*b+3*c*c+a*c

#Exercise

#Find the solution to the following problem:
#Minimize:  x + y + 2x^2 + 2xy + 2y^2
#x - y >= 3
#x + 2y >= 6
#3x - 4y = 10 
#4 <= x <= 20
#2 <= y <= 22