""" This module solves the solves the multi-bang control problem min 1/2 \|y-yd\|^2 + alpha/2 \int |u|^2 + beta \prod_{i=1}^d |u-u_i|_0 s.t. -\Delta y = u, u_1 <= u(x) <= u_d using the approach described in the paper "Multi-bang control of elliptic systems" by Christian Clason, and Karl Kunisch, see http://math.uni-graz.at/mobis/publications/SFB-Report-2013-006.pdf. Besides NumPy and SciPy, it requires the 'dolfin' module from the open source FEniCS Project (http://fenicsproject.org/). """ __author__ = "Christian Clason " __date__ = "March 27, 2013" import numpy as np import scipy.sparse as sp from scipy.sparse.linalg import spsolve from dolfin import * # problem parameters N = 128 # number of nodes per dimension maxit = 10 # max number of Newton steps alpha = 5e-3 # penalty parameter (L2) beta = 1e-3 # penalty parameter (L0) ui = np.linspace(-2,2,5) # vector of control states # construct scipy matrix from bilinear form a, boundary condition bc parameters['linear_algebra_backend'] = 'uBLAS' def assemble_csr(a,bc): row,col,val = assemble(a,bcs=bc).data() return sp.csr_matrix((val,col,row)) # setup grid, target, differential equation mesh = UnitSquareMesh(N,N) V = FunctionSpace(mesh,'CG',1) u = TrialFunction(V) v = TestFunction(V) # target yd = Expression('3/10.*pow(1-(6*x[0]-3),2)*exp(-pow((6*x[0]-3),2)\ - pow((6*x[1]-3)+1,2)) - (1./5*(6*x[0]-3) - pow((6*x[0]-3),3)\ - pow((6*x[1]-3),5))*exp(-pow((6*x[0]-3),2)-pow((6*x[1]-3),2))\ - 1./30*exp(-pow((6*x[1]-3)+1,2) - pow((6*x[1]-3),2))') plot(yd,mesh=mesh,title='target') # Laplace equation with homogeneous Dirichlet conditions a = dot(grad(u),grad(v))*dx bc = DirichletBC(V, Constant(0.0), lambda x,on_boundary: on_boundary) A = assemble_csr(a,bc) # stiffness matrix AT = assemble_csr(adjoint(a),bc) # stiffness matrix (adjoint) M = assemble_csr(u*v*dx,bc) # mass matrix z = assemble(yd*v*dx,bcs=bc) # load vector n = V.dim() # number of nodes d = ui.shape[0] # number of control states # initialization and preallocation y,p = np.zeros(n), np.zeros(n) # state, adjoint state He, DHe = np.zeros(n), np.zeros(n) # H_eps(p^k), D_N H_eps(p^k) As = np.zeros([n,4*d-2]) # active sets As_old = np.zeros([n,4*d-2]) sqab = sqrt(2*alpha*beta) # continuation: start with sufficiently small eps^0 eps = min(sqab,np.min(alpha/4*(ui[1:]-ui[:d-1])))*0.9 while eps > 1e-12: print '\nCompute solution for epsilon = %1.3e' % (eps) for it in xrange(maxit): # update active sets As_old[:] = As # P_i^eps As[:,0] = (p <= (alpha*ui[0] + sqab-eps)) & \ (p <= (alpha/2*(ui[0]+ui[1])-eps)) for i in xrange(1,d-1): As[:,i] = (abs(p-alpha*ui[i]) <= (sqab-eps)) & \ (p <= (alpha/2*(ui[i]+ui[i+1])-eps)) & \ (p >= (alpha/2*(ui[i-1]+ui[i])+eps)) As[:,d-1] = (p >= (alpha*ui[d-1] - sqab+eps)) & \ (p >= (alpha/2*(ui[d-2]+ui[d-1])+eps)) He[:] = np.dot(As[:,:d],ui) # P_i,i+1^eps for i in xrange(d-1): ind = d+1+i if (alpha/2*(ui[i+1]-ui[i]) <= sqab): As[:,ind] = abs(p-alpha/2*(ui[i]+ui[i+1])) < eps He += As[:,ind]/(2*eps)*(ui[i]*(alpha/2*(ui[i]+ui[i+1])+eps-p)+ ui[i+1]*(p-alpha/2*(ui[i]+ui[i+1])+eps)) DHe[:] = np.dot(As[:,d+1:2*d],np.diff(ui))/(2*eps) # P_i,0-^eps, P_i,0+^eps for i in xrange(1,d-1): ind = 2*d + 2*i-2 if (alpha/2*(ui[i]-ui[i-1]) > sqab): As[:,ind] = abs(p-alpha*ui[i]+sqab) < eps He += ((alpha*ui[i]-sqab-eps)*(alpha*ui[i]-sqab+eps-p)/alpha + ui[i]*(p-alpha*ui[i]+sqab+eps))*As[:,ind]/(2*eps) if (alpha/2*(ui[i+1]-ui[i]) > sqab): As[:,ind+1] = abs(p-alpha*ui[i]-sqab) < eps He += As[:,ind+1]/(2*eps)*(ui[i]*(alpha*ui[i]+sqab+eps-p) +\ (alpha*ui[i]+sqab+eps)*(p-alpha*ui[i]-sqab+eps)/alpha) DHe += (sqab+eps)/(2*eps*alpha)*np.sum(As[:,2*d:],1) # P_0^eps: all remaining points As[:,d] = 1-np.sum(As[:,:d],1)-np.sum(As[:,d+1:],1) He += As[:,d]*p/alpha DHe += As[:,d]/alpha # setup Newton matrix, right hand side IS = sp.spdiags(DHe,0,n,n,format='csr') H = sp.bmat([[M,AT],[A,-M*IS]],format='csr') g = np.hstack([z-M*y-AT*p,-A*y+M*He]) nr = np.linalg.norm(g) update = np.count_nonzero(As-As_old) print 'It# %d: update = %i,\t residual = %1.3e' % (it+1,update,nr) if update == 0: break # semismooth Newton step dy,dp = np.split(spsolve(H,g),2) y += dy p += dp # check convergence if it < maxit-1: # converged: accept iterate uvec = He regnodes = np.count_nonzero(As[:,d+1:]) print 'Solution has %i node(s) in regularized active sets' % (regnodes) if regnodes == 0: break # solution optimal: terminate else: eps/=10. # reduce eps, continue Newton else: # not converged: reject, terminate print 'Iterate rejected, returning u_eps for eps = %1.3e' %(eps*10) break # plot control ufun = Function(V) ufun.vector()[:] = uvec plot(ufun,title='control') interactive()