library(quadprog)
# Simulating data from a linear model
set.seed(2019)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 0.3 * x1 - 0.3 * x2 + rnorm(n, sd=0.5)
X <- cbind(x1,x2)
# estimating a linear model with no constraints
r <- lm(y ~ -1 + x1 + x2)

## adding equality constraints

# e.g. the coefficients should sum up to b
A <- matrix(1,2,1) # constraint matrix
b <- 0.2
# using solve.QP
eqR <- solve.QP(t(X)%*%X, t(X)%*%y, A, b, meq=1)
# manual implementation
LHS <- rbind(cbind(t(X)%*%X, A), cbind(t(A), 0))
RHS <- c(t(y)%*%X, b)
eqM <- solve(LHS, RHS)
# compare outcomes
cbind(eqR$solution, eqM[1:2])
         [,1]       [,2]
x1  0.4450044  0.4450044
x2 -0.2450044 -0.2450044
# check constraint
t(A)%*%eqR$solution; t(A)%*%eqM[1:2]
     [,1]
[1,]  0.2
     [,1]
[1,]  0.2

## adding inequality constraints

# e.g. the coefficients should be larger than b
A <- diag(2)
b <- c(0.452,-0.33)
ineqR <- solve.QP(t(X)%*%X, t(X)%*%y, A, b)
# check constraints
t(A)%*%ineqR$solution >= b
     [,1]
[1,] TRUE
[2,] TRUE
## manual
# here I was not able to reproduce manually ineqR
# meaning solely using solve()