最適配置の数理　朝倉書店　岡部篤行・鈴木敦夫先生　

```r

#p.57

#長方形領域(x:x1<=x<=x2,y:y1<=y<=y2)で人口密度がf(x,y)である場合の最適配置の位置を計算するアルゴリズム
#2重積分をモンテカルロ法により計算し、最急降下法で最適化

f=function(x,y){
return((1/2)*dnorm(x,2,1)*dnorm(y,2,1)+(1/2)*dnorm(x,4,1)*dnorm(y,4,1))  
}

x1=1
x2=5
y1=1
y2=6

#モンテカルロ積分によるメッシュ
N=10^3
xs=runif(N,x1,x2)
ys=runif(N,y1,y2)

t=function(xvec,yvec){
value=0
for(i in 1:length(yvec)){
for(j in 1:length(xvec)){
value=value+sum(sqrt((xvec[j]-xs)^2+(yvec[i]-ys)^2)*f(xs,ys))
}}    
return(value/(N))
}

#初期値
xpoint=runif(1,x1,x2)
ypoint=runif(1,y1,y2)

ite=3000
eta=10^(-1)
h=10^(-3)

for(l in 1:ite){

for(i in 1:length(xpoint)){
xpoint_sub=xpoint
xpoint_sub[i]=xpoint_sub[i]+h
xpoint[i]=xpoint[i]-eta*(t(xpoint_sub,ypoint)-t(xpoint,ypoint))/h

ypoint_sub=ypoint
ypoint_sub[i]=ypoint_sub[i]+h
ypoint[i]=ypoint[i]-eta*(t(xpoint,ypoint_sub)-t(xpoint,ypoint))/h

}  
    
print(t(xpoint,ypoint))  
  
}

```

```r

#p.98

#長方形領域(x:x11<=x<=x12,y:y11<=y<=y12)で人口密度がf1(x,y),長方形領域(x:x21<=x<=x22,y:y21<=y<=y22)で人口密度がf2(x,y)である場合の最適配置の位置を計算するアルゴリズム(2点の距離の制約条件付き)

f1=function(x,y){
return((1/2)*dnorm(x,2,1)*dnorm(y,2,1)+(1/2)*dnorm(x,4,1)*dnorm(y,4,1))  
}

f2=function(x,y){
return(dnorm(x,7,1)*dnorm(y,8,1))  
}

#f1(x,y)での領域
x11=1
x12=5
y11=1
y12=6

#f2(x,y)での領域
x21=5
x22=10
y21=1
y22=6

#モンテカルロ積分によるメッシュ
N=10^3
xs1=runif(N,x11,x12)
ys1=runif(N,y11,y12)
xs2=runif(N,x21,x22)
ys2=runif(N,y21,y22)

lam=1 #平滑化パラメータ(制約条件の強さ)
K=4.5 #2点の距離の制約条件

t=function(xvec,yvec){
value1=sum(sqrt((xvec[1]-xs1)^2+(yvec[1]-ys1)^2)*f1(xs1,ys1))
value2=sum(sqrt((xvec[2]-xs2)^2+(yvec[2]-ys2)^2)*f2(xs2,ys2))
value=(value1+value2)/N
distance=sqrt((xvec[1]-xvec[2])^2+(yvec[1]-yvec[2])^2)
if(distance>K){
value=value+lam*(distance-K)^2  
}
return(value)
}

#初期値
xpoint1=runif(1,x11,x12)
ypoint1=runif(1,y11,y12)
xpoint2=runif(1,x21,x22)
ypoint2=runif(1,y21,y22)
xpoint=c(xpoint1,xpoint2)
ypoint=c(ypoint1,ypoint2)


ite=10^4
eta=5*10^(-1)
h=10^(-3)

for(l in 1:ite){

for(i in 1:length(xpoint)){
#x座標の更新
xpoint_sub=xpoint
xpoint_sub[i]=xpoint_sub[i]+h
xpoint[i]=xpoint[i]-eta*(t(xpoint_sub,ypoint)-t(xpoint,ypoint))/h
#y座標の更新
ypoint_sub=ypoint
ypoint_sub[i]=ypoint_sub[i]+h
ypoint[i]=ypoint[i]-eta*(t(xpoint,ypoint_sub)-t(xpoint,ypoint))/h
}  
#コスト関数の値    
print(t(xpoint,ypoint))  
}

#計算後の結果
print(paste0("x点:",xpoint))
print(paste0("y点:",ypoint))
print(paste0("2点の距離:",sqrt((xpoint[1]-xpoint[2])^2+(ypoint[1]-ypoint[2])^2)))


```

```r

#最適配置を2つの領域で計算し、人口密度の最も過密な経路を探索するアルゴリズム

#長方形領域(x:x1<=x<=x2,y:y1<=y<=y2)で人口密度がf(x,y)である場合の最適配置の位置を計算するアルゴリズム
#2重積分をモンテカルロ法により計算し、最急降下法で最適化

#f=function(x,y){
#return((1/2)*dnorm(x,2,1)*dnorm(y,2,1)+(1/2)*dnorm(x,4,1)*dnorm(y,4,1))  
#}

f=function(x,y){
return(dnorm(x,7,1)*dnorm(y,8,1))  
}

#x1=1
#x2=5
#y1=1
#y2=6

x1=5
x2=10
y1=1
y2=6

#モンテカルロ積分によるメッシュ
N=10^3
xs=runif(N,x1,x2)
ys=runif(N,y1,y2)

t=function(xvec,yvec){
value=0
for(i in 1:length(yvec)){
for(j in 1:length(xvec)){
value=value+sum(sqrt((xvec[j]-xs)^2+(yvec[i]-ys)^2)*f(xs,ys))
}}    
return(value/(N))
}

#初期値
xpoint=runif(1,x1,x2)
ypoint=runif(1,y1,y2)

ite=3000
eta=10^(-1)
h=10^(-3)

for(l in 1:ite){

for(i in 1:length(xpoint)){
xpoint_sub=xpoint
xpoint_sub[i]=xpoint_sub[i]+h
xpoint[i]=xpoint[i]-eta*(t(xpoint_sub,ypoint)-t(xpoint,ypoint))/h

ypoint_sub=ypoint
ypoint_sub[i]=ypoint_sub[i]+h
ypoint[i]=ypoint[i]-eta*(t(xpoint,ypoint_sub)-t(xpoint,ypoint))/h

}  
    
print(t(xpoint,ypoint))  
  
}

#f1(x,y)
#> xpoint
#[1] 3.052578
#> ypoint
#[1] 3.07312

#f2(x,y)
#> xpoint
#[1] 7.988524
#> ypoint
#[1] 1.66433


#p.98(応用)

#長方形領域(x:x11<=x<=x12,y:y11<=y<=y12)で人口密度がf1(x,y),長方形領域(x:x21<=x<=x22,y:y21<=y<=y22)で人口密度がf2(x,y)である場合の2点の経路で人口密度がもっとも多いと思われる経路計算のアルゴリズム

f1=function(x,y){
return((1/2)*dnorm(x,2,1)*dnorm(y,2,1)+(1/2)*dnorm(x,4,1)*dnorm(y,4,1))  
}

f2=function(x,y){
return(dnorm(x,7,1)*dnorm(y,8,1))  
}

#f1(x,y)での領域
x11=1
x12=5
y11=1
y12=6

#f2(x,y)での領域
x21=5
x22=10
y21=1
y22=6

#モンテカルロ積分によるメッシュ
N=10^3
xs1=runif(N,x11,x12)
ys1=runif(N,y11,y12)
xs2=runif(N,x21,x22)
ys2=runif(N,y21,y22)

#初期値
xpoint1=3.052578 #f1(x,y)における最適な配置
ypoint1=3.07312 #f1(x,y)における最適な配置
xpoint2=7.988524 #f2(x,y)における最適な配置
ypoint2=1.66433 #f2(x,y)における最適な配置
dist=sqrt((xpoint1-xpoint2)^2+(ypoint1-ypoint2)^2)


mesh=0.1
Kmesh=seq(mesh,dist,mesh)
lam=1 #平滑化パラメータ(制約条件の強さ)

xpoints1=c(xpoint1)
ypoints1=c(ypoint1)
xpoints2=c(xpoint2)
ypoints2=c(ypoint2)

for(p in 1:length(Kmesh)){

K=Kmesh[p] #2点の距離の制約条件

t=function(xvec,yvec){
value1=sum(sqrt((xvec[1]-xs1)^2+(yvec[1]-ys1)^2)*f1(xs1,ys1))
value2=sum(sqrt((xvec[2]-xs2)^2+(yvec[2]-ys2)^2)*f2(xs2,ys2))
value=(value1+value2)/N
distance=sqrt((xvec[1]-xvec[2])^2+(yvec[1]-yvec[2])^2)
#if(distance>K){
value=value+lam*(distance-K)^2  
#}
return(value)
}

#初期値
xpoint=c(xpoint1,xpoint2)
ypoint=c(ypoint1,ypoint2)

ite=10^3
eta=10^(-1)
h=10^(-3)

for(l in 1:ite){

for(i in 1:length(xpoint)){
#x座標の更新
xpoint_sub=xpoint
xpoint_sub[i]=xpoint_sub[i]+h
xpoint[i]=xpoint[i]-eta*(t(xpoint_sub,ypoint)-t(xpoint,ypoint))/h
#y座標の更新
ypoint_sub=ypoint
ypoint_sub[i]=ypoint_sub[i]+h
ypoint[i]=ypoint[i]-eta*(t(xpoint,ypoint_sub)-t(xpoint,ypoint))/h
}  
#コスト関数の値    
#print(t(xpoint,ypoint))  
}

#計算後の結果
print(paste0("x点:",xpoint))
print(paste0("y点:",ypoint))
print(paste0("2点の距離:",sqrt((xpoint[1]-xpoint[2])^2+(ypoint[1]-ypoint[2])^2)))

xpoints1=c(xpoints1,xpoint[1])
ypoints1=c(ypoints1,ypoint[1])
xpoints2=c(xpoints2,xpoint[2])
ypoints2=c(ypoints2,ypoint[2])

}

xpoints=c(xpoints1,xpoints2)
ypoints=c(ypoints1,ypoints2)

#人口密度が最も多いと思われる2点の経路
plot(xpoints,ypoints)

```