restart; with(plots): with(LinearAlgebra): #exampleNo := 1: # Beispiel 1 - Ornstein-Uhlenbeck # Beispiel 2 - Geometric Brownian # Beispiel 3 - Shiryaev ### !!!!! ### ### Die folgende Zeile bei erstmaligem Ausfuehren auskommentieren und auch nur dann! Ansonsten muss das Worksheet ggf. neu geoeffnet werden. #currentdir("../../"); ### !!!!! ### ### In dem unterem Verzeichnis ggf. "ex3-shiryaev" in "ex1-ou" oder "ex2-gb" anpassen, je nach dem, welcher Prozess behandelt werden soll. ### Danach den darauffolgenden Ordner anpassen, je nach dem, welche numerische Simulation geplottet werden soll. currentdir("ut_N1_T005"); #currentdir(); readParams := proc() global exampleNo, xLeft, lengthX, mu_gamma, sigma, stateDim, controlDim, lambdaValue, alphaValue, mpcHorizon, mpcSteps, mpcStepsize, TT: local params: params := ImportVector("config.txt", skiplines=1): exampleNo := params[1]: xLeft := params[2]: lengthX := params[3]: mu_gamma := params[4]: sigma := params[5]: stateDim := convert(params[6],rational)+2: # +2 Randpunkte. controlDim := convert(params[7],rational): lambdaValue := convert(convert(params[10],rational),float,1): alphaValue := convert(params[11],rational): mpcHorizon := convert(params[14],rational): mpcSteps := convert(params[15],rational): mpcStepsize := convert(convert(params[16],rational),float,1): TT := mpcSteps * mpcStepsize: # Endzeit = MPC Sampling Rate * MPC Steps end proc: numericSol := proc() local coords, values, numericPlot; coords := Vector([seq(xLeft + ii*lengthX/(stateDim - 1), ii=0..(stateDim - 1))]); values := ImportMatrix("trajectory.txt"); print(values); if(exampleNo = 1) then numericPlot := display([pointplot(coords, Transpose(Row(values,1)), connect, legend="y_0", linestyle="dot", title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=0.0), titlefont=[Times,Bold,12]), seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y", title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=2..mpcSteps+1)], insequence=true, color=blue); elif(exampleNo = 2) then numericPlot := display([pointplot(coords, Transpose(Row(values,1)), connect, legend="y_0", linestyle="dot", title=typeset("Geometric Brownian\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=0.0), titlefont=[Times,Bold,12]), seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y", title=typeset("Geometric Brownian\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue,t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=2..mpcSteps+1)], insequence=true, color=blue); elif(exampleNo = 3) then numericPlot := display([pointplot(coords, Transpose(Row(values,1)), connect, legend="y_0", linestyle="dot", title=typeset("Shiryaev\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=0.0), titlefont=[Times,Bold,12]), seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y", title=typeset("Shiryaev\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=2..mpcSteps+1)], insequence=true, color=blue); fi: return numericPlot; end proc: adjointPlot := proc() local coords, values, adjointplot; coords := Vector([seq(xLeft + ii*lengthX/(stateDim - 1), ii=0..(stateDim - 1))]); values := ImportMatrix("adjTrajectory.txt"); print(values); if(exampleNo = 1) then adjointplot := display([seq(pointplot(coords, Transpose(Row(values,i)), connect, labels=[x,p(x)], title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=1..mpcSteps),plot(0, x=xLeft..xLeft+lengthX, labels=[x,p(x)], title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=TT), titlefont=[Times,Bold,12], color=green)], insequence=true, color=green); elif(exampleNo = 2) then adjointplot := display([seq(pointplot(coords, Transpose(Row(values,i)), connect, labels=[x,p(x)], title=typeset("Geometric Brownian\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=1..mpcSteps), plot(0, x=xLeft..xLeft+lengthX, labels=[x,p(x)], title=typeset("Geometric Brownian\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=TT), titlefont=[Times,Bold,12], color=green)], insequence=true, color=green); elif(exampleNo = 3) then adjointplot := display([seq(pointplot(coords, Transpose(Row(values,i)), connect, labels=[x,p(x)], title=typeset("Shiryaev\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=1..mpcSteps),plot(0, x=xLeft..xLeft+lengthX, labels=[x,p(x)], title=typeset("Shiryaev\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=TT), titlefont=[Times,Bold,12], color=green)], insequence=true, color=green); fi: return adjointplot; end proc: controlPlot := proc() local controls, controlTimes, controlPoints, controlplot, coords; if(controlDim = 1) then controls := Column(ImportMatrix("control.txt"),1): controlTimes := Vector([seq(floor(ii/2)*mpcStepsize, ii=1..2*mpcSteps)]); controlPoints := Vector([seq(controls[ceil(ii/2)], ii=1..2*mpcSteps)]); print(controls); if(exampleNo = 1) then controlplot := pointplot(controlTimes, controlPoints, connect, labels=[t,u(t)], title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12]); elif(exampleNo = 2) then controlplot := pointplot(controlTimes, controlPoints, connect, labels=[t,u(t)], title=typeset("Geometric Brownian\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12]); elif(exampleNo = 3) then controlplot := pointplot(controlTimes, controlPoints, connect, labels=[t,u(t)], title=typeset("Shiryaev\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12]); fi; elif(controlDim = stateDim-2) then coords := Vector([seq(xLeft + ii*lengthX/(stateDim - 1), ii=1..(stateDim - 2))]); controls := ImportMatrix("control.txt"); if(ColumnDimension(controls) > stateDim) then controls := controls[1..RowDimension(controls),1..ColumnDimension(controls)]; fi: #print(coords); print(controls); if(exampleNo = 1) then controlplot := display([seq(pointplot(coords, Transpose(Row(controls,i)), connect, labels=[x,u(x)], title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=1..mpcSteps), plot(0, x=xLeft..xLeft+lengthX, labels=[x,u(x)], title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=TT), titlefont=[Times,Bold,12], color=black)], insequence=true, color=black, view=[xLeft..xLeft+lengthX,default]); elif(exampleNo = 2) then controlplot := display([seq(pointplot(coords, Transpose(Row(controls,i)), connect, labels=[x,u(x)], title=typeset("Geometric Brownian\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=1..mpcSteps), plot(0, x=xLeft..xLeft+lengthX, labels=[x,u(x)], title=typeset("Geometric Brownian\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=TT), titlefont=[Times,Bold,12], color=black)], insequence=true, color=black, view=[xLeft..xLeft+lengthX,default]); elif(exampleNo = 3) then controlplot := display([seq(pointplot(coords, Transpose(Row(controls,i)), connect, labels=[x,u(x)], title=typeset("Shiryaev\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=(i-1)*mpcStepsize), titlefont=[Times,Bold,12]), i=1..mpcSteps), plot(0, x=xLeft..xLeft+lengthX, labels=[x,u(x)], title=typeset("Shiryaev\134n %1, %2, %3, %4, %5", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue, t=TT), titlefont=[Times,Bold,12], color=black)], insequence=true, color=black, view=[xLeft..xLeft+lengthX,default]); fi; fi; return controlplot; end proc: targetFctPlot := proc() local coords, values, targetfctplot; coords := Vector([seq(xLeft + ii*lengthX/(stateDim - 1), ii=0..(stateDim - 1))]); values := ImportMatrix("desiredTrajectory.txt"); print(values); if(exampleNo = 1) then targetfctplot := display([plot(0, x=xLeft..xLeft+lengthX, title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12], color=red),seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y_d", title=typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12], color=red, linestyle=dash), i=1..mpcSteps)], insequence=true, color=red); elif(exampleNo = 2) then targetfctplot := display([plot(0, x=xLeft..xLeft+lengthX, title=typeset("Geometric Brownian\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12], color=red), seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y_d", title=typeset("Geometric Brownian\134n %1, %2, %3 %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12], color=red, linestyle=dash), i=1..mpcSteps)], insequence=true, color=red); elif(exampleNo = 3) then targetfctplot := display([plot(0, x=xLeft..xLeft+lengthX, title=typeset("Shiryaev\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12], color=red), seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y_d", title=typeset("Shiryaev\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12], color=red, linestyle=dash), i=1..mpcSteps)], insequence=true, color=red); fi: return targetfctplot; end proc: numericSol3d := proc() local coords, values, numericPlot; coords := Vector([seq(xLeft + ii*lengthX/(stateDim - 1), ii=0..(stateDim - 1))]); values := ImportMatrix("trajectory.txt"); print(values); numericPlot := display([seq(pointplot(coords, Transpose(Row(values,i)), connect, legend="y"), i=1..mpcSteps+1)], insequence=true, color=blue); return numericPlot; end proc: readParams(); #mpcSteps := 143: totalCosts := ImportMatrix("totalcosts.txt", skiplines=1)[1,5]; exampleNo; xLeft; lengthX; mu_gamma; sigma; stateDim; controlDim; lambdaValue; alphaValue; mpcHorizon; mpcSteps; mpcStepsize; TT; #targetFct1 := (x,t) -> 1/(sqrt(2*Pi*0.2^2)) * exp(-(x-2*sin(Pi*t/5))^2 / (2*0.2^2)); #targetFct2 := (x,t) -> 1/(x*sqrt(2*Pi*0.1^2)) * exp(-(log(x) - (1 + sin(Pi*t/5)))^2 / (2*0.1^2)); #if(exampleNo = 1) then # #bild1 := display([seq(plot(targetFct1(x,i * mpcStepsize), x=xLeft..xLeft+lengthX, view=[default,0..4], color=red, linestyle="dash", legend="y_d", #title=typeset("Ornstein-Uhlenbeck\134n %1, %2, costs = %3", N=mpcHorizon, lambda=lambdaValue, totalCosts), titlefont=[Times,Bold,12]), i=0..mpcSteps)], #insequence=true): # bild1 := display([seq(plot(targetFct1(x,i * mpcStepsize), x=xLeft..xLeft+lengthX, view=[default,0..2], color=red, linestyle="dash", legend="y_d", title=#typeset("Ornstein-Uhlenbeck\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12]), i=0..mpcSteps)], insequence=#true): #elif(exampleNo = 2) then # bild1 := display([seq(plot(targetFct2(x,i * mpcStepsize), x=xLeft..xLeft+lengthX, view=[default,0..2], legend="y_d", color=red, linestyle="dash", title=#typeset("Geometric-Brownian\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12]), i=0..mpcSteps)], insequence=#true): #else # bild1 := display([seq(plot(targetFct2(x,i * mpcStepsize), x=xLeft..xLeft+lengthX, view=[default,0..2], legend="y_d", color=red, linestyle="dash", title=#typeset("Shiryaev\134n %1, %2, %3, %4", T=mpcStepsize, N=mpcHorizon, lambda=lambdaValue, alpha=alphaValue), titlefont=[Times,Bold,12]), i=0..mpcSteps)], insequence=true): #fi: bild1 := targetFctPlot(): bild2 := numericSol(): display(bild2, insequence=true); display(bild1, bild2); controlPlot(); adjointPlot();