function mineffort
% This function solves the regularized minimum effort problem
%   min 1/2 |y-z|^2 + alpha/2 c^2 + (beta c)/2 |u|^2 s.t. Ay=cu, |u(x)|<=1
% using semi-smooth Newton method and continuation in beta, with stopping
% criterion via model function (described in the paper
%      "A minimum effort optimal control problem for elliptic PDEs"      
% by Christian Clason, Kazufumi Ito, and Karl Kunisch, see
% http://math.uni-graz.at/mobis/publications/SFB-Report-2011-002.pdf.
%
% January 27, 2011          Christian Clason <christian.clason@uni-graz.at>
%                                            http://www.uni-graz.at/~clason

%% setup
N     = 128;        % grid size
alpha = 1e-3;       % control cost penalty (L^infty)
maxit = 10;         % max iterations in semismooth Newton method
bred  = 1/10;       % reduction factor for beta
eff   = 0.99;       % desired efficiency of stopping criterion

% coefficients of differential operator
nu = .1;            % diffusion
b  = -[0,1];        % convection

% setup grid on [-1,1]^2
[x,y] = meshgrid(linspace(-1,1,N));
h     = 2/N;
N2    = N^2;
alpha = alpha*N2;   % scale penalty to correct for finite dimensional norm

tplot = @(n,f,s) tplot_(n,f,s,N,x,y);

% setup differential operators
Dxx = gallery('poisson',N)*1/h^2;   % Laplace operator (negative)
e   = ones(N,1);
Dx  = 1/h*spdiags([-e e],-1:0,N,N); % derivative operator in 1D
Ix  = speye(N,N);

A   = nu*Dxx+b(1)*kron(Ix,Dx)+b(2)*kron(Dx,Ix);

% setup target
z = 1.*(((x-.5).^2+(y+.2).^2) < 1/32)+0.5*(((x+.2).^2+(y-.3).^2) < 1/16);

% precompute data, operators
Az = A*z(:); z = z(:);
AA = A*A';

%% solve

% initialize iterates, parameters
c = 0;              % optimal bound on u
p = zeros(N2,1);    % dual variable (Lagrange multiplier)

beta  = 1;          % regularization parameter
m0    = 1e99;       % predicted value of c_beta |u_beta|^2 at beta = 0
unorm = 0;          % norm of computed control u_beta

% continuation in beta
while (c*unorm < eff*m0) && (beta > 1e-12)
    
    fprintf('\nCompute minimizer for beta=%1.3e:\n',beta);
    
    cold = c;
    pold = p;
    
    unold = unorm;
    
    % semismooth Newton method
    iter   = 0;
    update = 1e99;
    
    while (update) && (iter < maxit)
        
        % active sets
        ap = p >  beta;
        am = p < -beta;
        as = ap+am;
        is = 1-as;       % inactive set
        
        sbp = ap - am + is.*(p/beta);                  % sign_beta(p)
        nbp = as'*(abs(p)-beta/2)+is'*(p.^2)/(2*beta); % |p|_beta
        
        % Hessian
        Lcc = alpha;
        Lcp = -sbp';
        Lpc = sbp;
        Lpp = AA + spdiags(c/beta*is,0,N2,N2);
        
        H = [Lcc Lcp; Lpc Lpp];
        
        % gradient
        Lc  = alpha*c - nbp;
        Lp  = AA*p + c*sbp - Az;
        rhs = -[Lc;Lp];
        nr  = norm(rhs(:));      % norm of residual
        
        % Newton step
        dx = H\rhs;
        
        c = c+dx(1);
        p = p+dx(2:end);
        
        % terminate if active sets coincide
        if iter == 0
            update = 1e99;
        else
            update = nnz((ap-ap_old) | (am-am_old));
            fprintf('It# %i: update = %i,\t residual = %1.3e\n',iter,update,nr);
        end
        ap_old = ap;
        am_old = am;
        
        iter = iter+1;
    end
    
    % steer continuation in beta
    if update                   % semi-smooth Newton iteration did not converge
        c    = cold;            % reject step
        p    = pold;
        bred = (bred)^(1/2);    % reduce step size for beta
        fprintf('Reduction factor increased to %1.3e\n',bred);
        beta = bold*bred;       % try again with larger beta
    else                        % accept step
        u = ((ap-am)+(is.*p/beta));
        y = A\(c*u);
        Jval  = (1/2 * norm(y-z)^2 + alpha/2*c^2)/N2;
        unorm = norm(u)^2/N^2;
        
        fprintf('\nc_beta = %g,\t J(y_beta,c_beta) = %g\n',c, Jval);
        
        % build model function for termination of continuation
        if beta < 1
            % interpolation points: last two iterates
            b1 = bold;
            c1 = cold*unold;
            b2 = beta;
            c2 = c*unorm;
            
            k1 = c2*c1*(-2*b1*c2*b2+c1*b1^2+2*b1*(c2*c1*(b1-b2)^2)^(1/2)+...
                c2*b2^2-2*b2*b1*c1-2*b2*(c2*c1*(b1-b2)^2)^(1/2)+...
                c2*b1^2+b2^2*c1)/(-c1+c2)^2;
            k2 = (-c2*b2+b1*c1+(c2*c1*(b1-b2)^2)^(1/2))/(-c1+c2);
            
            mfun = @(x) k1./(x+k2).^2;    % model function
            m0   = mfun(0);               % intercept
        end
        
        bold = beta;
        beta = beta*bred;                 % reduce beta
    end
    
end

%% plot
tplot(1,z,'target');
tplot(2,y,'optimal state y^*')
tplot(3,p,'optimal Lagrange multiplier p^*')
tplot(4,c*u,'optimal control (c^*u^*)')

end % main function

function tplot_(n,f,s,N,x,y)
figure(n); surf(x,y,reshape(f,N,N));
shading interp; camlight headlight; lighting gouraud; alpha(0.8);
title(s); xlabel('x'); ylabel('y');
set(gca,'FontSize',12); drawnow;
end % tplot_ function