### user equilibrium
 # リンクパフォーマンス関数の宣言
pt1<-function(x1) 10*(1+0.15*(x1/2)^4)
pt2<-function(x2) 20*(1+0.15*(x2/4)^4)
pt3<-function(x3) 25*(1+0.15*(x3/3)^4)


#リンクフロー数の初期値宣言
x1<-0
x2<-0
x3<-0
X<-10

#リンクフロー数に基づくリンク旅行時間の計算結果表示→どの経路の旅行時間が最も小さいか
t1<-pt1(x1=x1)
t2<-pt2(x2=x2)
t3<-pt3(x3=x3)

#リンクフローに基づくリンクフロー数の割り当て(all or nothing)→最も短いところに全ての交通量を流す
x1<-(t1==min(t1,t2,t3))*(X)
x2<-(t2==min(t1,t2,t3))*(X)
x3<-(t3==min(t1,t2,t3))*(X)
#リンクフロー数結果表示
x1
x2
x3
#リンク旅行時間更新→旅行時間に基づき配分した結果旅行時間が変化するので再度計算
t1<-pt1(x1=x1)
t2<-pt2(x2=x2)
t3<-pt3(x3=x3)
#リンクフロー数に基づくリンク旅行時間の計算結果表示
t1
t2
t3

#収束判定値の初期化
d <- Inf  #収束判定値の初期値は無限大に設定
while(d > 0.1){       #dが0.0001よりも大きい限り{}の操作を継続
#補助リンクフロー数→とりあえず流した結果最も短い旅行時間のところへ新たに交通量を配分する準備
y1<-(t1==min(t1,t2,t3))*(X)
y2<-(t2==min(t1,t2,t3))*(X)
y3<-(t3==min(t1,t2,t3))*(X)

#補助リンクフロー数表示
y1
y2
y3

#ラインサーチ→どれだけの交通量を配分すればいいか
library(cubature)  #計算に必要なライブラリーの呼び出し(先に積分計算をしたものを計算した方がいいけど)
fl<-function(alpha) adaptIntegrate(pt1,0,x1+alpha*(y1-x1))$integral+adaptIntegrate(pt2,0,x2+alpha*(y2-x2))$integral+adaptIntegrate(pt3,0,x3+alpha*(y3-x3))$integral

ans<-optimize(fl,lower=0,upper=1)
alpha<-ans$minimum
alpha
#リンクフロー数の更新
x1<-x1+alpha*(y1-x1)
x2<-x2+alpha*(y2-x2)
x3<-x3+alpha*(y3-x3)

#リンクフロー数表示
x1
x2
x3

#リンク旅行時間更新
t1<-pt1(x1=x1)
t2<-pt2(x2=x2)
t3<-pt3(x3=x3)

t1
t2
t3
d1<-min(t1,t2,t3)
#補助リンクフロー数
y1<-(t1==min(t1,t2,t3))*(X)
y2<-(t2==min(t1,t2,t3))*(X)
y3<-(t3==min(t1,t2,t3))*(X)

#補助リンクフロー数表示
y1
y2
y3
 #ラインサーチ
fl<-function(alpha) adaptIntegrate(pt1,0,x1+alpha*(y1-x1))$integral+adaptIntegrate(pt2,0,x2+alpha*(y2-x2))$integral+adaptIntegrate(pt3,0,x3+alpha*(y3-x3))$integral

ans<-optimize(fl,lower=0,upper=1)
alpha<-ans$minimum
alpha
#リンクフロー数の更新
x1<-x1+alpha*(y1-x1)
x2<-x2+alpha*(y2-x2)
x3<-x3+alpha*(y3-x3)
 #リンクフロー数表示
x1
x2
x3

#リンク旅行時間更新
t1<-pt1(x1=x1)
t2<-pt2(x2=x2)
t3<-pt3(x3=x3)

t1
t2
t3


#目的関数値
ans$objective
d2<-min(t1,t2,t3)
d<-(abs(d1-d2)/d2)
}
#結果の出力
x1
x2
x3
t1
t2
t3
d