0

我正在尝试为一些控制问题编写马尔可夫链近似。但是我在 R 中有以下错误，我在 Stackoverflow 中检查了类似的问题，但仍然不知道如何解决它。任何帮助将不胜感激。

该错误来自我想在 for 循环中找到所有“u”中的最小值的地方。

具体来说，在uit-for-loop, 对于每一个 nextuit我可以得到一个新的单个值（我想）temp，并想将它与单个值变量存储的临时最小值进行比较vmin。这就是if-else句子中的想法。

最好跳过参数设置和初始化过程。

#----- parameters ------
xleft=0; xright=10
yleft=0; yright=10
h=0.01
Nx=(xright-xleft)/h
Ns=2
Nu=11; hu=0.2
la=0.1
qMainDiag=c(-0.5,-0.5)
qSubDiag=c(0.5,0.5)
alpha=c(0.2,0.25)
beta=c(0.35,0.2)
a=c(0.6,0.8)
b=c(0.5,0.3)
c=c(0.45,0.5)
d=c(0.65,0.8)
tol=10^(-8)
maxitr=10000
#---- Initialization -----
Vold=array(0,dim=c(Nx+1,Nx+1,Ns))
Vnew=array(0,dim=c(Nx+1,Nx+1,Ns))
Uopt=array(0,dim=c(Nx+1,Nx+1,Ns))
for(r in 1:Ns){
  for(i in 1:(Nx+1)){
    for(j in 1:(Nx+1)){
      Vold[i,j,r]=1
    }
  }
}
#---- iteration ----
for(n in 1:maxitr){
  for(r in 1:Ns){
    # inner of O
    for(i in 2:Nx){
      for(j in 2:Nx){
        vInt=0
        for(it in 1:(min(i,j)+1)){
          vInt=vInt+Vold[i-it+1,j-it+1,r]*0.1*exp(-0.1*(it-1)*h)*h
        }
        # For each u, want to find the minimum temp value and its u.
        for(uit in 1:Nu){
          x=xleft+(i-1)*h; y=yleft+(j-1)*h
          u=hu*(uit-1)
          Xi11=(alpha[r]*x)^2; Xi22=(beta[r]*y)^2
          f1=x*(a[r]-b[r]*y+u); f2=y*(-c[r]+d[r]*x+u)
          g=1+r*(x+y)*(1+u^2)
          Qh=(Xi11+Xi22)+h*(abs(f1)+abs(f2))+h-(h^2)*qMainDiag[r]
          dlt=(h*h)/Qh
          pforward=0.5*(Xi11+2*h*max(f1,0.0))/Qh
          pback=0.5*(Xi11+2*h*max(-f1,0.0))/Qh
          pup=0.5*(Xi22+2*h*max(f2,0.0))/Qh
          pdown=0.5*(Xi22+2*h*max(-f2,0.0))/Qh
          pswitch=(h*h*qSubDiag[r])/Qh
          pstay=h/Qh
          temp=(1-la*dlt)*(pforward*Vold[i+1,j,r]+pback*Vold[i-1,j,r]
                           +pup*Vold[i,j+1,r]+pdown*Vold[i,j-1,r]
                           +pswitch*Vold[i,j,3-r]
                           +pstay*Vold[i,j,r])+la*dlt*vInt+dlt*g
          # find the minimal value (Here is the spot!!!)
          if(uit==1){
            vmin=temp; umin=u
          }else if(temp<vmin){
            vmin=temp; umin=u
          }
        }
        Vnew[i,j,r]=vmin
        Uopt[i,j,r]=umin
      }
    }
    errormax=max(abs(Vold-Vnew))
    print(n)
    print(errormax)
    Vold=Vnew
    if(errormax<tol){
      break
    }
  }
}
4

0 回答 0