Fusion-Power-Plant-Framework / bluemira

Bluemira is an integrated inter-disciplinary design tool for future fusion reactors. It incorporates several modules, some of which rely on other codes, to carry out a range of typical conceptual fusion reactor design activities.
https://bluemira.readthedocs.io/
GNU Lesser General Public License v2.1
67 stars 16 forks source link

Current Midplance Constraint identifies eq as PlasmaCoil object, instead of an Equilibrium one #3077

Closed athoynilimanew closed 6 months ago

athoynilimanew commented 8 months ago

Describe the bug

I have been trying to use the CurrentMidplanceConstraint function from bluemira.equilibria.optimisation.constraint_funcs, to use with a TikhonovCurrentCOP

But I am getting the following error:

File [~/WorkDir/code/bluemira/bluemira/equilibria/optimisation/constraint_funcs.py:247], in CurrentMidplanceConstraint.f_constraint(self, vector)
    [245] """Constraint function"""
    [246] self.eq.coilset.get_control_coils().current = vector * self.scale
--> [247] lcfs = self.eq.get_LCFS()
    [248] if self.inboard:
    [249]     return self.radius - min(lcfs.x)

AttributeError: 'PlasmaCoil' object has no attribute 'get_LCFS'

Steps to reproduce

This is how I built a constraint class inherited from UpdateableConstraint

class CurrentMidplanceConstraintClass(UpdateableConstraint):
    """
    Constraint to constrain the inboard or outboard midplane
    of the plasma during optimisation.

    Can be used as a standalone constraint for use in an optimisation problem.
    """

    def __init__(
        self,
        coilset: CoilSet,
        radius: float,
        scale: float,
        inboard: bool,
        tolerance: Union[float, np.ndarray] = 1.0e-6,
        constraint_type: str = "inequality",
    ):  
        n_coils = coilset.n_coils()
        if is_num(tolerance):
            tolerance = tolerance * np.ones(n_coils)
        if len(tolerance) != n_coils:
            raise ValueError(
                "Tolerance vector length not equal to the number of coils."
            )

        # arguments needed for CurrentMidplanceConstraint function
        self._args = {
            "radius": radius,
            "scale": scale,
            "inboard": inboard,
        }
        self.f_constraint_type = constraint_type
        self.tolerance = tolerance

    def prepare(self, equilibrium: Equilibrium, I_not_dI=False, fixed_coils=True):
        """
        Prepare the constraint for use in an equilibrium optimisation problem.

        This is called in each iteration. Better to put parameters
        that needed to be calculated dynamically during optimisation here

        """
        if I_not_dI:
            equilibrium = _get_dummy_equilibrium(equilibrium)

        self._args["eq"] = equilibrium

    def control_response(self, coilset: CoilSet):
        """
        Calculate control response of a CoilSet to the constraint.
        """
        pass

    def evaluate(self, equilibrium: Equilibrium):
        """
        Calculate the value of the constraint in an Equilibrium.
        """
        pass

    def f_constraint(self) -> CurrentMidplanceConstraint:
        """The numerical non-linear part of the constraint."""
        f_constraint = CurrentMidplanceConstraint(**self._args)
        f_constraint.constraint_type = self.f_constraint_type
        return f_constraint

and this is how I use it:

current_midplane_constraint = CurrentMidplanceConstraintClass(
    coilset=coilset,
    radius=1.0,
    scale=1.0,
    inboard=True,
    tolerance=1e-4
)
opt_problem = TikhonovCurrentCOP(
    coilset,
    eq,
    magnetic_targets,
    gamma=1e-8,
    max_currents=3.0e7,
    opt_algorithm="SLSQP",
    opt_conditions={"max_eval": 100},
    opt_parameters={},
    constraints=[current_midplane_constraint],
)
constrained_iterator = PicardIterator(
    eq,
    opt_problem,
    fixed_coils=True,
    plot=False,
    relaxation=0.3,
    convergence=DudsonConvergence(1e-4),
)

constrained_iterator()

given coilset and eq defined before.

Expected behaviour

I expected the constrained_iterator() to run. As defined in the bluemira.equilibria.optimisation.constraint_funcs, CurrentMidplanceConstraint should know that self.eq is an Equilibriumobject, not a PlasmaCoil.

Screenshots / Tracebacks

Environment (please complete the following information)

Additional context

Discussed this with @oliverfunk. We tried to modify CurrentMidplanceConstraint to have LCFSas input, instead of Equilibrium. The same error as above was observed.

You can see the old current midplane constraint definition in here: old coilset optimisation example

CoronelBuendia commented 8 months ago

Ah, I can explain this a little bit. When doing a current vector optimisation (as opposed to dcurrent vector optimisation - see _get_dummy_equilibrium) we trick the constraint so that it splits out the plasma background response from the controlled coil response.

The PlasmaCoil is masquerading deliberately as an Equilibrium. We'd never had any constraints that required methods on Equilibrium that were not also on PlasmaCoil.

The simple answer is that LCFS is an optional input and get called up when I_not_dI == True, because in the PlasmaCoil case the concept of an LCFS does not really make sense. I am not sure that will do what you need. Personally, in order to comply more easily with our constraint framework, I would use a normalised flux equaity constraint with a slightly looser than normal tolerance, but those units would sadly not be [m].

Bit of a side note: I have never really understood this constraint... It is similar to what one might call a "gap" constraint in normal parlance, which is something that I have never implemented. Technically, I think they are only valid about a "linearised" plasma state, for which I'm told you need to converge representative equilibria to get a good idea of the linearised control response of a gap constraints. I'm actually thinking about trying to do an approximate gap constraint response function, but never seem to find the time...

je-cook commented 6 months ago

As this is by design, closing