View source for Source Inversion (FEniCS/Casadi)

This page contains a discretized version of the mixed integer PDE constrained optimization problem [[Source Inversion]] in [http://www.casadi.org CasADi 2.4.3] format. The weak version of the problem was discretized using finite element methods (via [http://www.fenicsproject.org FEniCS]) using first-degree Lagrangian elements on an equidistant mesh with <math>128</math> cells.

=== CasADi ===

<source lang="Python">
from dolfin import *

import math
import numpy
import scipy, scipy.sparse
import casadi


# Simple Branch and Bound solver
def solve_bnb(solver, root, integers):
    # Initialize the node queue
    Q = [(numpy.inf, root)]

    # Global solution data
    x_opt = None
    f_opt = numpy.inf

    # Continuous solution data
    x0 = None
    f0 = numpy.inf

    # Main loop
    seqnum = 0
    while Q:
        seqnum += 1

        # Extract node with lowest lower bound and check if fathoming is necessary
        # (should not happen)
        node = Q.pop()
        if f_opt < numpy.inf and node[0] >= f_opt:
            print "Node %4d (%4d left): %9g >= %9g - PRE-FATHOM" % (seqnum, len(Q), node[0], f_opt)
            continue
        node = node[1]

        # Solve the node
        result = solver(node)

        # Local solution data
        x = result['x']
        f = result['f'].get()[0]

        # Store continuous solution data if none has been stored yet
        if x0 is None:
            x0 = casadi.DMatrix(x)
            f0 = f

        # Check if branch can be fathomed
        if f >= f_opt:
            print "Node %4d (%4d left): %9g >= %9g - POST-FATHOM" % (seqnum, len(Q), f, f_opt)
            continue

        # Check for violations of integrality (fixed tolerance 1e-5)
        viol = [abs(casadi.floor(x[i] + 0.5) - x[i]) for i in integers]
        idx  = [(integers[i], viol[i]) for i in range(len(integers)) if viol[i] > 1e-5]
        if not idx:
            # Register new global solution
            x_opt = x
            f_opt = f

            # Cull the branch and bound tree
            pre  = len(Q)
            Q    = [n for n in Q if n[0] < f_opt]
            post = len(Q)

            print "Node %4d (%4d left): f_opt = %9g - *** SOLUTION *** (%d culled)" % (seqnum, len(Q), f, pre - post)
        else:
            # Branch on first violation
            idx = idx[0][0]

            # Generate two new nodes (could reuse node structure)
            ln = {
                'x0':     casadi.DMatrix(x),
                'lbx':    node['lbx'],
                'ubx':    casadi.DMatrix(node['ubx']),
                'lbg':    node['lbg'],
                'ubg':    node['ubg'],
                'lam_x0': casadi.DMatrix(result['lam_x']),
                'lam_g0': casadi.DMatrix(result['lam_g'])
            }
            un = {
                'x0':     casadi.DMatrix(x),
                'lbx':    casadi.DMatrix(node['lbx']),
                'ubx':    node['ubx'],
                'lbg':    node['lbg'],
                'ubg':    node['ubg'],
                'lam_x0': casadi.DMatrix(result['lam_x']),
                'lam_g0': casadi.DMatrix(result['lam_g'])
            }

            ln['ubx'][idx] = casadi.floor(x[idx])
            un['lbx'][idx] = casadi.ceil(x[idx])

            lower_node = (f, ln)
            upper_node = (f, un)

            # Insert new nodes in queue (inefficient for large queues)
            Q.extend([lower_node, upper_node])
            Q.sort(cmp=lambda x,y: cmp(y[0], x[0]))
            print "Node %4d (%4d left): %9g - BRANCH ON %d" % (seqnum, len(Q), f, idx)

    return {
        'x0': x0,
        'f0': f0,
        'x': x_opt,
        'f': f_opt
    }


# Converts FEniCS matrix to scipy.sparse.coo_matrix
def matrix_to_coo(A, shape=None):
    nnz = A.nnz()
    col = numpy.empty(nnz, dtype='uint64')
    row = numpy.empty(nnz, dtype='uint64')
    val = numpy.empty(nnz)

    if shape is None:
        shape = (A.size(0), A.size(1))

    it = 0
    for i in range(A.size(0)):
        col_local, val_local = A.getrow(i)
        nnz_local = col_local.shape[0]

        col[it:it+nnz_local] = col_local[:]
        row[it:it+nnz_local] = i
        val[it:it+nnz_local] = val_local[:]
        it                  += nnz_local

    return scipy.sparse.coo_matrix((val, (row, col)), shape=shape)


# Returns FEniCS expression for elementary source term
def source_term(x0):
    return Expression("100.0 * exp(-pow(x[0] - x0, 2) / 0.02)", x0=x0)

# Grid resolutions
nx = 128
nu = 8

# Generate mesh and function space
mesh = UnitIntervalMesh(nx)
V    = FunctionSpace(mesh, 'CG', 1)

# Define the bilinear form and assemble it
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx + u * v * ds
A = assemble(a)

# Create the linear form for the exact source term and perform assembly
f = source_term(3.0 / math.pi) + source_term(0.5)
L = f * v * dx
b = assemble(L)

# Solve for the reference solution
u_ref = Function(V)
solve(A, u_ref.vector(), b)

# Define the control grid and the elementary source terms
ctrl_grid = numpy.linspace(0.0, 1.0, num=nu+1)[:-1]
ctrl_grid += 0.5 * ctrl_grid[1]
src_expr  = [source_term(x0) for x0 in ctrl_grid]

# Assemble the linear forms corresponding to the elementary source terms and
# construct the second half of the coefficient matrix
b = [assemble(-f * v * dx) for f in src_expr]
B = scipy.column_stack([v.array() for v in b])

# Construct the full coefficient matrix
C = scipy.sparse.vstack([
    scipy.sparse.hstack([matrix_to_coo(A), B]),
    scipy.sparse.coo_matrix(([1.0]*nu, ([0]*nu, [V.dim() + i for i in range(nu)])), shape=(1, V.dim() + nu))
])

# Transfer the optimization constraints to CasADi
C  = casadi.DMatrix(C)
lb = casadi.DMatrix.zeros(C.shape[0], 1)
ub = casadi.DMatrix.zeros(C.shape[0], 1)

lb[-1] = 0.0
ub[-1] = 5.0

lbx = numpy.array([-numpy.inf]*V.dim() + [0.0]*nu)
ubx = numpy.array([ numpy.inf]*V.dim() + [1.0]*nu)

# Define the objective functional and assemble it
u     = Function(V)
J     = (u_ref - u)**2 * dx
dJdu  = derivative(J, u)
dJdu2 = derivative(dJdu, u)

H = casadi.DMatrix(matrix_to_coo(assemble(dJdu2), shape=(V.dim() + nu, V.dim() + nu)))
g = casadi.DMatrix(numpy.concatenate((assemble(dJdu).array(), [0.0]*nu)))
c = assemble(J)

# Build the CasADi NLP
x   = casadi.MX.sym('X', V.dim() + nu)
nlp = casadi.MXFunction('nlp',
    casadi.nlpIn(x=x),
    casadi.nlpOut(
        f=0.5*casadi.quad_form(x, H) + casadi.inner_prod(g, x) + c,
        g=casadi.mul(C, x)
    ))

# Build NLP solver
solver = casadi.NlpSolver('solver', 'ipopt', nlp, {
    'print_level': 0,
    'print_time':  False
})
result = solve_bnb(solver, {
    'lbx': lbx,
    'ubx': ubx,
    'lbg': lb,
    'ubg': ub
}, range(V.dim(), V.dim() + nu))

# Evaluate result
print ""
print "----------------------------"
print " f = %g" % result['f']

u.vector()[:] = result['x'][:V.dim()].toArray()[:,0]
</source>

== Results ==

Node solutions are calculated using IPOPT (3.12, default settings, using linear solver MUMPS on Ubuntu 15.10 for x86-64 with Linux 4.2.0-19-generic, Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz, 16 GB RAM). The objective function value is <math>0.858837</math>. The integer solution is found after processing 25 nodes.

== Solution ==

The following is a Gnuplot compatible tabular version of the control data determined by the solver:

<pre>
#       x         w        w0
   0.0625         0 2.91549e-08
   0.1875         0 2.75108e-08
   0.3125         0 6.746e-08
   0.4375         1   0.46917
   0.5625         0  0.543817
   0.6875         0 1.33205e-07
   0.8125         0 1.51069e-07
   0.9375         1  0.897707
</pre>

The following is a Gnuplot compatible tabular version of the solution functions:

<pre>
#       x         u        u0       u_d
        0    19.845   18.7521   18.7182
0.0078125   20.0001   18.8986   18.8645
 0.015625   20.1551   19.0451   19.0107
0.0234375   20.3101   19.1916   19.1569
  0.03125   20.4652   19.3381   19.3032
0.0390625   20.6202   19.4846   19.4494
 0.046875   20.7752   19.6311   19.5956
0.0546875   20.9303   19.7776   19.7419
   0.0625   21.0853   19.9241   19.8881
0.0703125   21.2403   20.0706   20.0343
 0.078125   21.3953   20.2171   20.1806
0.0859375   21.5503   20.3635   20.3268
  0.09375   21.7053     20.51    20.473
 0.101562   21.8603   20.6565   20.6193
 0.109375   22.0153   20.8029   20.7655
 0.117188   22.1702   20.9494   20.9117
    0.125   22.3251   21.0958   21.0579
 0.132812   22.4799   21.2422   21.2042
 0.140625   22.6347   21.3886   21.3504
 0.148438   22.7894   21.5349   21.4966
  0.15625    22.944   21.6812   21.6428
 0.164062   23.0985   21.8275   21.7889
 0.171875   23.2528   21.9737   21.9351
 0.179688    23.407   22.1197   22.0812
   0.1875    23.561   22.2657   22.2273
 0.195312   23.7146   22.4116   22.3733
 0.203125    23.868   22.5573   22.5193
 0.210938   24.0209   22.7028   22.6652
  0.21875   24.1734    22.848    22.811
 0.226562   24.3254    22.993   22.9567
 0.234375   24.4766   23.1377   23.1022
 0.242188   24.6271    23.282   23.2476
     0.25   24.7767   23.4259   23.3927
 0.257812   24.9252   23.5692   23.5376
 0.265625   25.0725   23.7119   23.6822
 0.273438   25.2184    23.854   23.8263
  0.28125   25.3627   23.9952     23.97
 0.289062   25.5052   24.1355   24.1131
 0.296875   25.6457   24.2748   24.2556
 0.304688   25.7839    24.413   24.3972
   0.3125   25.9196   24.5498    24.538
 0.320312   26.0525   24.6852   24.6777
 0.328125   26.1823   24.8189   24.8162
 0.335938   26.3088   24.9509   24.9533
  0.34375   26.4316   25.0809   25.0889
 0.351562   26.5505   25.2087   25.2226
 0.359375   26.6651   25.3342   25.3542
 0.367188   26.7753   25.4572   25.4836
    0.375   26.8807   25.5774   25.6105
 0.382812   26.9811   25.6947   25.7346
 0.390625   27.0762   25.8089   25.8556
 0.398438   27.1659   25.9198   25.9733
  0.40625   27.2499   26.0271   26.0873
 0.414062   27.3281   26.1308   26.1973
 0.421875   27.4004   26.2305   26.3032
 0.429688   27.4666   26.3262   26.4046
   0.4375   27.5268   26.4177   26.5012
 0.445312   27.5809   26.5047   26.5927
 0.453125   27.6289   26.5873   26.6791
 0.460938   27.6708   26.6652   26.7599
  0.46875   27.7069   26.7384   26.8351
 0.476562   27.7371   26.8066   26.9045
 0.484375   27.7617     26.87    26.968
 0.492188   27.7808   26.9283   27.0255
      0.5   27.7947   26.9816   27.0768
 0.507812   27.8035   27.0298   27.1221
 0.515625   27.8076   27.0729   27.1613
 0.523438   27.8072   27.1109   27.1945
  0.53125   27.8026   27.1438   27.2217
 0.539062    27.794   27.1718   27.2431
 0.546875   27.7818   27.1948   27.2589
 0.554688   27.7662    27.213   27.2692
   0.5625   27.7476   27.2264   27.2743
 0.570312   27.7261   27.2352   27.2743
 0.578125   27.7021   27.2395   27.2696
 0.585938   27.6759   27.2394   27.2603
  0.59375   27.6476   27.2351   27.2469
 0.601562   27.6174   27.2269   27.2295
 0.609375   27.5857   27.2147   27.2084
 0.617188   27.5525    27.199    27.184
    0.625   27.5181   27.1797   27.1565
 0.632812   27.4826   27.1573   27.1261
 0.640625   27.4461   27.1317   27.0932
 0.648438   27.4087   27.1032    27.058
  0.65625   27.3706   27.0721   27.0207
 0.664062   27.3318   27.0385   26.9815
 0.671875   27.2924   27.0025   26.9406
 0.679688   27.2525   26.9644   26.8983
   0.6875   27.2119   26.9242   26.8545
 0.695312   27.1709   26.8821   26.8096
 0.703125   27.1293   26.8383   26.7635
 0.710938   27.0871   26.7928   26.7164
  0.71875   27.0443   26.7457   26.6683
 0.726562   27.0008   26.6971   26.6193
 0.734375   26.9566   26.6469   26.5693
 0.742188   26.9115   26.5953   26.5184
     0.75   26.8654   26.5422   26.4666
 0.757812   26.8183   26.4875   26.4137
 0.765625   26.7699   26.4312   26.3598
 0.773438     26.72   26.3732   26.3046
  0.78125   26.6686   26.3135   26.2481
 0.789062   26.6153   26.2517   26.1902
 0.796875     26.56    26.188   26.1306
 0.804688   26.5024   26.1219   26.0691
   0.8125   26.4423   26.0534   26.0057
 0.820312   26.3794   25.9822   25.9399
 0.828125   26.3134   25.9082   25.8717
 0.835938    26.244   25.8311   25.8007
  0.84375    26.171   25.7506   25.7267
 0.851562    26.094   25.6665   25.6494
 0.859375   26.0129   25.5785   25.5685
 0.867188   25.9272   25.4865   25.4837
    0.875   25.8368   25.3902   25.3947
 0.882812   25.7413   25.2893   25.3013
 0.890625   25.6406   25.1837   25.2032
 0.898438   25.5345   25.0732   25.1002
  0.90625   25.4227   24.9576   24.9919
 0.914062    25.305   24.8368   24.8782
 0.921875   25.1815   24.7107   24.7589
 0.929688   25.0519   24.5791   24.6338
   0.9375   24.9162   24.4421   24.5029
 0.945312   24.7745   24.2996   24.3659
 0.953125   24.6267   24.1516   24.2228
 0.960938   24.4728   23.9982   24.0737
  0.96875    24.313   23.8395   23.9184
 0.976562   24.1474   23.6756   23.7572
 0.984375   23.9762   23.5067   23.5899
 0.992188   23.7995   23.3328   23.4169
        1   23.6175   23.1542   23.2381
</pre>
 
[[Category:Casadi]]