ファジィ動的計画法　槇書店　小田中敏夫先生

```r

#p.7 ピカールの方法

#曲線y=1,x=0
b=1
a=0

#初期値
y=1
X=0.5

ite=100
yvec=c(y)

for(i in 1:ite){

#シンプソンの第一公式による積分値の計算(一変数関数の積分)

#関数の定義：
func=function(x){
z<-1+x*y
return(z)
}

#a:積分区間の始点、a+H:積分区間の終点、n:積分区間の分割数
a=0;H=X;n=1000

#重合公式の計算
Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));
for(j in 1:(n-2)){
  #奇数番目は4、偶数番目は2になる係数
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))
}

y=b+Y
yvec=c(yvec,y)
}

#収束する
print(yvec)

```

```r

#p.14 費用が2*(x-p)^2,廃棄が20円の費用の場合にnカ月間の費用最小とする生産量を求める

p=c(210,250,300,200,400)
s=c(180,195,220,210,300)

f=function(x,P,S){
return(2*(x-P)^2+20*(x-S))
}

fi=function(x){
fval=c() 
for(j in 1:length(x)){  
if(j>1){
value=0  
for(i in 2:j){
value=value+f(x[j],p[i],s[i])+f(x[j],p[i-1],s[i-1])
}}else{
value=f(x[1],p[1],s[1])  
}
fval=c(fval,value)
}
return(fval)  
}


ite=150
eta=0.01
X=p+10
h=0.001

for(l in 1:ite){
  
dX=c()  
for(j in 1:length(X)){
Xsub=X;Xsub[j]=Xsub[j]+h
dX=c(dX,(fi(Xsub)-fi(X))[j]/h)
}  

X=X-eta*dX
#売上予想分は需要となるため、生産量がそれを下回らないようにする。
X=ifelse(X<s,s,X)
print(fi(X))
  
}


```


```r

#p.31 1.3.5 折れ線近似
#制約条件が必要

g=function(x){
return(sin(x))  
}

CalcInt=function(u1,u2){
n=1500
A=(u2^3-u1^3)/3
B=(u2^2-u1^2)/2
C=u2-u1
X=u1+(u2-u1)*seq(1,n,1)/n
D=(u2-u1)*sum(g(X)^2)/n
E=(u2-u1)*sum(X*g(X)^2)/n
f=(u2-u1)*sum(g(X))/n  
return(c(A,B,C,D,E,f))
}

linearParam=function(vec){
A=vec[1];B=vec[2];C=vec[3];D=vec[4];E=vec[5];f=vec[6]
alpha=det(matrix(c(A,B,E,f),nrow=2))/det(matrix(c(A,B,B,C),nrow=2))
beta=det(matrix(c(E,f,B,C),nrow=2))/det(matrix(c(A,B,B,C),nrow=2))
return(c(alpha,beta))
}

IntCost=function(vec){
A=vec[1];B=vec[2];C=vec[3];D=vec[4];E=vec[5];f=vec[6]  
return(det(matrix(c(D,E,f,E,A,B,f,D,C),nrow=3))/det(matrix(c(A,B,B,C),nrow=2)))}

#全計算領域
a=1
b=pi

Cost=function(U){
value=0
for(i in 1:length(U)){
if(i==1){
param=CalcInt(a,a+U[1])
value=value+IntCost(param)}
if((i!=1)&(i!=length(U))){
param=CalcInt(a+sum(U[1:i]),a+sum(U[1:(i+1)]))
value=value+IntCost(param)}
if(i==length(U)){
param=CalcInt(a+sum(U[1:length(U)]),b)
value=value+IntCost(param)}
}
return(value)
}

newton=function(U,mue,lambda){
alpha=1
value=c()
for(i in 1:length(U)){
if(i==1){
value1=(g(a+U[1])-sum(linearParam(CalcInt(a,a+U[1]))*c(1,a+U[1])))^2 
value2=(g(a+U[1])-sum(linearParam(CalcInt(a+U[1],a+sum(U[1:2])))*c(1,a+U[1])))^2
value=c(value,value1-value2)}
if((i!=1)&(i!=length(U))){
value1=(g(a+sum(U[1:i]))-sum(linearParam(CalcInt(a+sum(U[1:(i-1)]),a+sum(U[1:i])))*c(1,a+sum(U[1:i]))))^2 
value2=(g(a+sum(U[1:i]))-sum(linearParam(CalcInt(a+sum(U[1:i]),a+sum(U[1:(i+1)])))*c(1,a+sum(U[1:i]))))^2
value=c(value,value1-value2)}
if(i==length(U)){
value1=(g(a+sum(U[1:i]))-sum(linearParam(CalcInt(a+sum(U[1:(i-1)]),a+sum(U[1:i])))*c(1,a+sum(U[1:i]))))^2 
value2=(g(a+sum(U[1:i]))-sum(linearParam(CalcInt(a+sum(U[1:i]),b))*c(1,a+sum(U[1:i]))))^2
value=c(value,value1-value2)  
}
}
value=value-mue/U+2*(sum(U)-b-a)
#value=c(value,-sum(log(U)),(sum(U)-(b-a))^2)
return(value) 
}

#newton-method
N=5
u=rep((b-a)/N,N-1)
z=rep(1,N-1)
mu=100
lambda=10^2
tau=0.999
ite=10^4
eta=10^(-4)
h=10^(-4)
const=10

while(const>10^(-4))({
ddL=array(0,dim=c(length(u),length(u)))
for(j in 1:length(u)){
usub=u
usub[j]=usub[j]+h
ddL[,j]=(newton(usub,mu,lambda)-newton(u,mu,lambda))/(h)
}  
dg=2*(sum(u)-(b-a))*rep(1,length(u))
dL=newton(u,mu,lambda)

mat1=cbind(ddL,t(t(c(dg))),-diag(1,length(dg)))
mat2=rbind(c(dg,rep(0,length(dg)+1)),cbind(diag(z),rep(0,length(dg)),diag(u)))
mat=rbind(mat1,mat2)
vec=c(-dL,(sum(u)-(b-a))^2,u*z-mu)
const=sqrt(sum(vec^2))*mu
duzy=solve(mat)%*%vec
u=u+eta*duzy[1:length(u)]
#u=(b-a)*u/(sum(u))
lambda=lambda+eta*duzy[length(u)+1]
z=z+eta*duzy[-c(1:(length(u)+1))]
mu=mu*tau

print(Cost(u))
})

```

```r

#最短経路問題　p.122

N=9
mat=matrix(c(0,9,3,rep(NA,6),9,0,2,4,rep(NA,5),3,2,0,8,2,7,rep(NA,3),NA,1,8,0,NA,NA,5,NA,NA,NA,NA,2,NA,0,NA,6,4,NA,NA,NA,7,NA,NA,0,NA,2,NA,rep(NA,3),5,6,NA,0,3,2,rep(NA,4),4,NA,3,0,6,rep(NA,6),2,6,0),nrow=n)
diag(mat)=0
mat[is.na(mat)>0]=10^4

#f[n,i]:iから出発してnステップで1に到達する最短距離
f=matrix(0,nrow=N,ncol=N)
f[1,]=mat[1,]
f[1,1]=0

for(n in 2:N) { 
	for(i in 1:N) {
	  #状態iに1から出発して次のステップで到達する時間を入力。
		f[n,i] = mat[i,1] + f[n-1,1] 
			for(j in 2:N){
			  #状態iに1以外から出発して次のステップで到達する最小時間を入力。
			  f[n,i] = min(f[n,i], mat[i,j]+f[n-1,j])
			}
	}
}

```