install.packages("Rsolnp") library(Rsolnp) #solnp function sig3<-c(2.5,3.0,5.0,16.0,26.06,36.25) sig1<-c(17.0,27.88,40.0,88.13,93.75,112.5) Paa<-101.25; Philimit<-0.01; Phiav<-29.29/180*pi PhiBeg<-Phiav-0.2; if (PhiBeg<0) PhiBeg<-0.0; PhiEnd<-Phiav+0.2 PhiBa<-seq(PhiBeg,PhiEnd,0.002) for (jPhi in 1:length(PhiBa)) { PhiCur<-PhiBa[jPhi] sigBa<-(sig1+sig3)/2- (sig1-sig3)/2*sin(PhiCur) tauBa<- (sig1-sig3)/2*cos(PhiCur) ObjFn <- function(par){ #1 A; 2 n; 3 T x1 <- par[1]; x2 <- par[2]; x3 <- par[3] tauCur<-Paa*x1*(sigBa/Paa+x3)^x2 sum(abs(tauCur-tauBa)^2)} x0 <- c(1,0.8,0.2) Res<-solnp(x0, fun = ObjFn,LB=c(0,0.5,0),UB=c(1000,1.0,500)) tanPhi<-Res$pars[1] *Res$pars[2] /((sigBa/Paa+Res$pars[3])^(1-Res$pars[2])) PhiNew<-atan(tanPhi); deltaPhi<-max(PhiNew-PhiCur) if (deltaPhi < Philimit) break } aBa<-Res$pars[1]; nBa<-Res$pars[2]; TBa<-Res$pars[3]; SOSba<-min(Res$values) SOSba