Merge pull request #17 from bendudson/optimise

Optimise equilibria

11-optimise-coils.py

@@ -0,0 +1,73 @@
+#!/usr/bin/env python
+
+import freegs
+
+#########################################
+# Create the machine, which specifies coil locations
+# and equilibrium, specifying the domain to solve over
+
+tokamak = freegs.machine.TestTokamak()
+
+eq = freegs.Equilibrium(tokamak=tokamak,
+                        Rmin=0.1, Rmax=2.0,    # Radial domain
+                        Zmin=-1.0, Zmax=1.0,   # Height range
+                        nx=65, ny=65,          # Number of grid points
+                        boundary=freegs.boundary.freeBoundaryHagenow)  # Boundary condition
+
+
+#########################################
+# Plasma profiles
+
+profiles = freegs.jtor.ConstrainPaxisIp(1e3, # Plasma pressure on axis [Pascals]
+                                        2e5, # Plasma current [Amps]
+                                        2.0) # Vacuum f=R*Bt
+
+#########################################
+# Coil current constraints
+#
+# Specify locations of the X-points
+# to use to constrain coil currents
+
+xpoints = [(1.1, -0.6),   # (R,Z) locations of X-points
+           (1.1, 0.6)]
+
+isoflux = [(1.1,-0.6, 1.1,0.6), # (R1,Z1, R2,Z2) pair of locations
+           (1.7, 0.0, 0.84, 0.0)]
+
+constrain = freegs.control.constrain(xpoints=xpoints, isoflux=isoflux, gamma = 1e-17)
+
+#########################################
+# Nonlinear solve
+
+freegs.solve(eq,          # The equilibrium to adjust
+             profiles,    # The toroidal current profile function
+             constrain)   # Constraint function to set coil currents
+
+
+
+# Currents in the coils
+tokamak.printCurrents()
+
+# Forces on the coils
+eq.printForces()
+
+############################
+# Optimise
+# Minimise the maximum force on the coils, while avoiding intersection of the LCFS and walls
+# by modifying the radius of the P2U and P2L coils.
+
+from freegs import optimise as opt
+
+best_eq = opt.optimise(eq,  # Starting equilibrium
+                       # List of controls
+                       [opt.CoilRadius("P2U"),
+                        opt.CoilRadius("P2L"), opt.CoilHeight("P2L")],
+                       # The function to minimise
+                       opt.weighted_sum(opt.max_coil_force, opt.no_wall_intersection),
+                       N=10,  # Number of solutions in each generation
+                       maxgen=20, # How many generations
+                       monitor=opt.PlotMonitor())  # Plot the best in each generation
+
+# Forces on the coils
+best_eq.printForces()
+best_eq.plot()

docs/index.rst

@@ -14,6 +14,7 @@ Welcome to FreeGS's documentation!
    input_and_output
    diagnostics
    tests
+   optimisation
 
 .. automodule:: freegs
 

docs/optimisation.rst

@@ -0,0 +1,102 @@
+.. _optimisation
+
+Optimisation
+============
+
+This is an experimental feature which is at an early stage of development.  The
+aim is to enable equilibria to be automatically optimised. This has the
+following components:
+
+#. Measures, a quantity which measures how "good" a solution is. Typically the
+   aim is to minimise this quantity, so I suppose it's really a measure of how
+   bad the solution is.
+#. Controls, quantities which can be changed. These could be machine parameters
+   such as coil locations, constraints like X-point location, or plasma profiles
+   such as poloidal beta or plasma current.
+#. An algorithm which modifies the controls and finds the best equilibrium
+   according to the measure it's given. At the moment the method used is
+   Differential Evolution.
+
+
+Differential Evolution
+----------------------
+
+`Differential Evolution <https://en.wikipedia.org/wiki/Differential_evolution>`_ is a type
+of stochastic search, similar to Genetic Algorithms, generally well suited to problems
+involving continuously varying parameters.
+
+The implementation of the algorithm is in ``freegs.optimiser``. It is generic,
+in that it operates on objects but does not need to know any details of what
+those objects are. To modify objects a list of ``controls`` are passed to the
+optimiser, each of which can set and get a value.  To score each object a
+``measure`` function is needed, which takes an object as input and returns a
+value. The optimiser works to minimise this value.
+
+An example which uses the optimisation method is in the ``freegs`` directory.
+This optimises a quadratic in 2D rather than tokamak equilibria. 100 generations
+are run, with 10 solutions (sometimes called agents) in each generation.  Run
+this example with the command:
+
+::
+
+   python test_optimiser.py
+
+This should produce the figure below. The red point is the best solution at each
+generation; black points are the other points in that generation. Faded colors
+(light red, grey) are used to show previous generations. It can be seen that the
+points are clustered around the starting solution, as the agents spread out, and
+then around the solution as the agents converge to the minimum.
+
+.. image:: optimiser.gif
+   :alt: Optimisation of a quadratic. The minimum is at (1,2), and starting point is at (0.1,0.1). Points mark the solutions tested at all generations.
+
+Optimising tokamak equilibria
+-----------------------------
+
+The code specific to FreeGS optimisation is in
+``freegs.optimise``. This includes controls which modify aspects of
+the tokamak or equilibrium, and measures which can be combined to
+specify the quantities which should be optimised.
+
++-------------------------------+----------------------------------------------------+
+|   Control                     |    Description                                     |
++===============================+====================================================+
+| CoilRadius(name [, min, max]) | Modify coil radius, given name and optional limits |
++-------------------------------+----------------------------------------------------+
+| CoilHeight(name [, min, max]) | Modify coil height                                 |
++-------------------------------+----------------------------------------------------+
+
+
++-------------------------------+----------------------------------------------------+
+|   Measure function            |   Description                                      |
++===============================+====================================================+
+| max_abs_coil_current          | Maximum current in any coil circuit                |
++-------------------------------+----------------------------------------------------+
+| max_coil_force                | Modify coil height                                 |
++-------------------------------+----------------------------------------------------+
+| no_wall_intersection          | Prevent intersections of wall and LCFS             |
++-------------------------------+----------------------------------------------------+
+
+Each measure function takes an Equilibrium as input, and returns a
+value. These can be combined in a weighted sum using
+``optimise.weighted_sum``, or by passing your own function to
+``optimise``.
+
+The example ``11-optimise-coils.py`` uses the following code to reduce
+the maximum force on the coils, while avoiding wall intersections::
+
+  from freegs import optimise as opt
+
+  best_eq = opt.optimise(eq,  # Starting equilibrium
+                         # List of controls
+                         [opt.CoilRadius("P2U"),
+                         opt.CoilRadius("P2L"), opt.CoilHeight("P2L")],
+                         # The function to minimise
+                         opt.weighted_sum(opt.max_coil_force, opt.no_wall_intersection),
+                         N=10,  # Number of solutions in each generation
+                         maxgen=20, # How many generations
+                         monitor=opt.PlotMonitor())  # Plot the best in each generation
+
+The monitor should be a callable (here it is an object of class ``PlotMonitor``), which
+is called after each generation. This is used to update a plot showing
+the best solution in each generation, and save the figure to file.

docs/tests.rst

@@ -1,6 +1,17 @@
 Tests
 =====
 
+Unit tests
+----------
+
+Unit testing is done with `pytest <https://docs.pytest.org/en/latest/>`_. Run the tests
+in the ``freegs`` directory with:
+
+::
+
+   pytest
+   
+
 Convergence test
 ----------------
 

freegs/README.md

@@ -3,6 +3,13 @@ FreeGS library
 
 Tests and examples using components of the library
 
+Testing
+-------
+
+To run the unit tests, install [pytest](https://docs.pytest.org/en/latest/) and then run:
+
+    $ pytest
+
 Multigrid solver
 ----------------
 

freegs/control.py

@@ -118,9 +118,12 @@ def __call__(self, eq):
         # Calculate the change in coil current
         current_change =  dot( inv(dot(transpose(A), A) + self.gamma**2 * eye(ncontrols)), 
                                dot(transpose(A),b))
-        print("Current changes: " + str(current_change))
+        #print("Current changes: " + str(current_change))
         tokamak.controlAdjust(current_change)
 
+        # Ensure that the last constraint used is set in the Equilibrium
+        eq._constraints = self
+
     def plot(self, axis=None, show=True):
         """
         Plots constraints used for coil current control
@@ -205,6 +208,9 @@ def __call__(self, eq):
         end_currents, _ = optimize.leastsq(self.psi_difference, start_currents, args=(eq,))
 
         tokamak.setControlCurrents(end_currents)
+
+        # Ensure that the last constraint used is set in the Equilibrium
+        eq._constraints = self
 
     def psi_difference(self, currents, eq):
         """
@@ -248,6 +254,9 @@ def __call__(self, eq):
 
         tokamak.setControlCurrents(end_currents)
 
+        # Ensure that the last constraint used is set in the Equilibrium
+        eq._constraints = self
+
     def psinorm_difference(self, currents, eq):
         """
         Difference between normalised psi from equilibrium with the given currents

freegs/equilibrium.py

@@ -41,7 +41,9 @@ class Equilibrium:
 
     _applyBoundary()
     _solver - Grad-Shafranov elliptic solver
-
+    _profiles     An object which calculates the toroidal current
+    _constraints  Control system which adjusts coil currents to meet constraints
+                  e.g. X-point location and flux values
     """
 
     def __init__(self, tokamak=machine.EmptyTokamak(),

freegs/jtor.py

@@ -241,7 +241,7 @@ def pshape(psinorm):
         L = self.Ip/I_R - LBeta0*(IR/I_R - 1)
         Beta0 = LBeta0 / L
 
-        print("Constraints: L = %e, Beta0 = %e" % (L, Beta0))
+        #print("Constraints: L = %e, Beta0 = %e" % (L, Beta0))
 
         # Toroidal current
         Jtor = L * (Beta0*R/self.Raxis + (1-Beta0)*self.Raxis/R)*jtorshape
@@ -372,7 +372,7 @@ def Jtor(self, R, Z, psi, psi_bndry=None):
         L = self.Ip/I_R - LBeta0*(IR/I_R - 1)
         Beta0 = LBeta0 / L
 
-        print("Constraints: L = %e, Beta0 = %e" % (L, Beta0))
+        #print("Constraints: L = %e, Beta0 = %e" % (L, Beta0))
 
         # Toroidal current
         Jtor = L * (Beta0*R/self.Raxis + (1-Beta0)*self.Raxis/R)*jtorshape

freegs/machine.py

@@ -736,7 +736,15 @@ def getForces(self, equilibrium = None):
         for label, coil in self.coils:
             forces[label] = coil.getForces(equilibrium)
         return forces
-
+
+    def getCurrents(self):
+        """
+        Returns a dictionary of coil label -> current in Amps
+        """
+        currents = {}
+        for label, coil in self.coils:
+            currents[label] = coil.current
+        return currents
 
 def EmptyTokamak():
     """