function oc_obstacle %% CODE FOR OPTIMAL CONTROL PROBLEM EXAMPLES IN PAPER BY ITO & KUNISCH %% % Authour : Serbiniyaz Anyyeva % Email: selbiniyaz(at)yahoo.com clear all; clc; beta=0.1; c=0.1; N=2^6; h=1/N; [X1,X2] = meshgrid(0:h:1, 0:h:1); %% Creating Discrete Laplacian for 2Dim G=zeros(N+1, N+1); numbers_v=1:1:(N-1)*(N-1); G(2:N,2:N)=reshape(numbers_v, N-1,N-1); LAPL_2 = delsq(G)/h/h; % -Delta (- laplacian) %% Meshing M=(N-1)*(N-1); [X1_int,X2_int] = meshgrid(h:h:1-h, h:h:1-h); % coordinates of the interior grid points Yd=reshape(yd(X1_int,X2_int), M,1); % vertical vector of yd grid values F=reshape(f(X1_int,X2_int), M,1); % vertical vector of f grid values PSI=reshape(psi(X1_int,X2_int), M,1); % vertical vector of psi grid values clear X1_int X2_int % Id=speye(M,M); Ze=sparse(zeros(M,M)); %% % % Solving unconstrained elliptic optimal control problem for initialization % $$\min 1/2*\|y-y_d\|^2+\beta/2*u^2 subj. \Delta y=u+f;$$ % A = horzcat(vertcat(Id,LAPL_2),vertcat(LAPL_2,-Id/beta)); V = [Yd; F]; SOL=A\V; clear A_init V_init % Y=SOL(1:M); P=SOL(M+1:2*M); % lambda = max(zeros(M,1),-LAPL_2*Y + P/beta+F); lambda_grid = G; S = max(zeros(M,1),lambda+c*(Y-PSI)) > 0; I=not(S); % %% Active set strategy loop_it=true; iter=0; while loop_it iter=iter+1; A=horzcat(vertcat(LAPL_2(I,I),Id(I,I)),vertcat(-Id(I,I)/beta,LAPL_2(I,I))); n_S=sum(S==1,1); n_I=M-n_S; P(S) = 0; Y(S)=PSI(S); V=vertcat(F(I)-LAPL_2(I,S)*Y(S),Yd(I)); SOL=A\V; Y(I)=SOL(1:n_I); P(I)=SOL(n_I+1:2*n_I); lambda(I)=0; lambda(S)=F(S)-(LAPL_2(S,S)*Y(S)+LAPL_2(S,I)*Y(I)); S1= max(zeros(M,1),lambda+c*(Y-PSI)) > 0; loop_it = not(and(all(S==S1),iter>1)); if (loop_it==false) disp('STOP condition is satisfied') end S=S1; I=not(S); if(iter==500) disp('Maximal number of iterations reached; breaking...') break end end %% Postprocessing Y_grid = G; Y_grid(G>0) = full(Y); figure(1); subplot(2,2,1); surf(X1,X2,Y_grid); title ('Y ') ; clear Y_grid; % U_grid = G; U_grid(G>0) = full(P/beta); subplot(2,2,2); surf(X1,X2,U_grid); title ('U ') ; clear U_grid; % lambda_grid = G; lambda_grid(G>0) = full(lambda); subplot(2,2,3); surf(X1,X2,lambda_grid); title ('\lambda '); clear lambda_grid; return function fun = f( X,Y ) fun=0.1*ones(size(X)); % fun = X-1/2; % fun=10^(-10)*ones(size(x)); return function fun = yd( X,Y ) % fun=zeros(size(X)); fun=5*X+Y-1; % fun=0.5*ones(size(X)); return function fun = psi( X,Y ) % fun=0*X; fun = 4*(X.*(X-1)+Y.*(Y-1))+1.5; % fun=ones(size(X)); return