最优化问题是普遍存在的,以前上运筹学课的时候也接触过最优化相关的问题,当时主要是理论课,并且关注的重点是单纯形法、运输问题以及图论等,这里指的最优化是指函数的最优化,即函数的极值,由于寻找一个局部最优比寻找全局最优要简单得多,所以这里的最优解也是指的局部最优解。
- 牛顿最优化方法
仅给出代码,公式什么的。。。我不知道博客园怎么插入公式,,,,
newton <- function(f3, x0, tol=1e-9, n.max = 100){
x <- x0
f3.x <- f3(x)
n <- 0
while((abs(f3.x[2]) > tol) & (n < n.max)){
x <- x - f3.x[2] / f3.x[3]
f3.x <- f3(x)
n <- n + 1
}
if(n == n.max){
cat("newton failed to converge\n")
}
else{
return(x)
}
} gamma.2.3 <- function(x){
if(x < 0) return(c(0,0,0))
if(x == 0) return(c(0,0,NaN))
y <- exp(-2*x)
return(c(4*x^2*y, 8*x*(1-x)*y, 8*(1-2*x^2)*y))
} x0 <- seq(0,10,0.01)
x <- c()
for(i in 1:length(x0)){
x[i] <- newton(gamma.2.3,x0[i])
}
代码在各种初始值下取到的极值点,在某些点可能会出现明显错误的极值点,因此使用的时候应该谨慎,多试试几组值。
- 黄金分割法
黄金分割法只能在一维情况下使用,但是不必知道函数的导数。其核心思想与夹逼求根的方法是一样的,即不断缩小区间范围。
为了加速运算,取了含1+ρ的项。
gsection <- function(ftn, x.l, x.r, x.m, tol=1e-9){
###黄金分割率
gr1 <- 1 + (1+sqrt(5))/2
f.l <- ftn(x.l)
f.r <- ftn(x.r)
f.m <- ftn(x.m)
while((x.r - x.l) > tol){
if((x.r - x.m) > (x.m - x.l)){
y <- x.m + (x.r - x.m)/grl
f.y <- ftn(y)
if(f.y >= f.m){
x.l <- x.m
f.l <- f.m
x.m <- y
f.m <- f.y
}
else{
x.r <- y
f.r <- f.y
}
}else{
y <- x.m - (x.m - x.l)/grl
f.y <- ftn(y)
if(f.y >= f.m){
x.r <- x.m
f.r <- f.m
x.m <- y
f.m <- f.y
}
else{
x.l <- y
f.l <- f.y
}
}
}
return(x.m)
}
f.4 <- function(x){return((x-1.7)^2+1)}
gsection(f.4, 1.5, 2, 1.8)