# 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 |p(x_)-y_i|.
# Let us constraining by min r, with r>=Au-y and r>=y-Au.
u0=c(3,0,1,0,0,0,0,100)  #p(x)=x^2+3 is an initial guess.
m=length(x)
n=length(u0)-1

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

Id

B=matrix(0,2*m,n+1)
B[1:m,1:n]=-A
B[(m+1):(2*m),1:n]=A
B[,n+1]=Id
B

yy=Id
yy[1:m]=-y
yy[(m+1):(2*m)]=y
## Solves linear and quadratic programming problems
## but needs a feasible starting value
#
# from example(solve.QP) in 'quadprog'
# 
fQP <- function(u) {u[length(u)]}
Amat       <- B
bvec       <- yy

gQP <- function(u) {v=0*(1:length(u))
    v[length(u)]=1
    v
}
####
result=constrOptim(u0, fQP,gQP, ui = B, 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)