# Test points 
x=c(1,2,3,4,5,6,7)
y=c(2,3,-1,3,2,-1,1)

plot(x,y)

# If p(x)=a+bx+x^2+dx^3+ex^4+fx^5+gx^6
# Let A=[x_i^(j-1)] and u=(a,b,...,f,g). We need to minimize (Au-y)^*(Au-y).
# Let us constraining by p(x)>=f(x) which means, in this case that Au>=y.
u0=c(3,0,1,0,0,0,0)  #p(x)=x^2+3 is an initial guess.
m=length(x)
n=length(u0)

A=matrix(0,m,n) 
for (i in 1:m){
for (j in 1:n){A[i,j]=x[i]^(j-1)}
}
A

## Solves linear and quadratic programming problems
## but needs a feasible starting value
#
# from example(solve.QP) in 'quadprog'
# 
fQP <- function(u) {t(A%*%u-y)%*%(A%*%u-y)}
Amat       <- A
bvec       <- y
# derivative
gQP <- function(u) {2*t(A)%*%(A%*%u-y)}
result=constrOptim(u0, fQP, gQP, ui = A, ci = bvec)
result
ur=result$par

poli<-function(t){
p=ur[1]
for (i in 1:(n-1)){p=p+ur[i+1]*t^i}
return(p)
}

t=seq(1,7,by=0.1);
z=c(poli(t[1]))
for ( i in 2:length(t)){z[i]=c(poli(t[i]))}
plot(t,z,'l')
points(x,y)