<jupyter_start><jupyter_text>Tutorial 2: Betz LimitNow that we have ran a simple calculator model using WISDEM, let's look at OpenMDAO. [OpenMDAO](http://openmdao.org) is the code that connects the various components of turbine models into a cohesive whole that can be optimized in systems engineering problems. WISDEM uses OpenMDAO to build up modular *components* and *groups* of components to represent a wind turbine.  Fortunately, OpenMDAO already provides some excellenet training examples on their [website](http://openmdao.org).  This tutorial is based on the OpenMDAO example, [Optimizing an Actuator Disk Model to Find Betz Limit for Wind Turbines](http://openmdao.org/twodocs/versions/latest/examples/betz_limit/betz.html), which we have extracted and added some additional commentary. The aim of this tutorial is to summarize the key points you'll use to create basic WISDEM models.  For those interested in WISDEM development, getting comfortable with all of the core OpenMDAO training examples is strongly encouraged.A classic problem of wind energy engineering is the Betz limit. This is the theoretical upper limit of power that can be extracted from wind by an idealized rotor. While a full explanation is beyond the scope of this tutorial, it is explained in many places online and in textbooks. One such explanation is on Wikipedia [https://en.wikipedia.org/wiki/Betz%27s_law](https://en.wikipedia.org/wiki/Betz%27s_law) . Problem formulationAccording to the Betz limit, the maximum power a turbine can extract from wind is:$$ C_p = \frac{16}{27} \approx 0.593 $$Where \\(C_p\\) is the coefficient of power.We will compute this limit using OpenMDAO by optimizing the power coefficient as a function of the induction factor (ratio of rotor plane velocity to freestream velocity) and a model of an idealized rotor using an actuator disk.Here is our actuator disc:![actuator disc](img/actuator_disc.png)Where \\(V_u\\) freestream air velocity, upstream of rotor, \\(V_r\\) is air velocity at rotor exit plane and \\(V_d\\) is slipstream air velocity downstream of rotor, all measured in \\(\frac{m}{s}\\).There are few other variables we'll have:- \\(a\\): Induced Velocity Factor- *Area*: Rotor disc area in \\(m^2\\)- *thrust*: Thrust produced by the rotor in N- \\(C_t\\): Thrust coefficient- *power*: Power produced by rotor in *W*- \\(\rho\\): Air density in \\(kg /m^3\\)Before we start in on the source code, let's look at a few key snippets of the code OpenMDAO implementationFirst we need to import OpenMDAO<jupyter_code>import openmdao.api as om<jupyter_output><empty_output><jupyter_text>Now we can make an `ActuatorDisc` class that models the actuator disc for the optimization.<jupyter_code>class ActuatorDisc(om.ExplicitComponent):
    def setup(self):
        # Inputs into the the model
        self.add_input('a', 0.5, desc='Indcued velocity factor')
        self.add_input('Area', 10.0, units='m**2', desc='Rotor disc area')
        self.add_input('rho', 1.225, units='kg/m**3', desc='Air density')
        self.add_input('Vu', 10.0, units='m/s', desc='Freestream air velocity, upstream of rotor')

        # Outputs
        self.add_output('Vr', 0.0, units='m/s', desc='Air velocity at rotor exit plane')
        self.add_output('Vd', 0.0, units='m/s', desc='Slipstream air velocity, downstream of rotor')
        self.add_output('Ct', 0.0, desc='Thrust coefficient')
        self.add_output('Cp', 0.0, desc='Power coefficient')
        self.add_output('power', 0.0, units='W', desc='Power produced by the rotor')
        self.add_output('thrust', 0.0, units='m/s')

        self.declare_partials('Vr', ['a', 'Vu'])
        self.declare_partials('Vd', 'a')
        self.declare_partials('Ct', 'a')
        self.declare_partials('thrust', ['a', 'Area', 'rho', 'Vu'])
        self.declare_partials('Cp', 'a')
        self.declare_partials('power', ['a', 'Area', 'rho', 'Vu'])

    def compute(self, inputs, outputs):
        a = inputs['a']
        Vu = inputs['Vu']
        rho = inputs['rho']
        Area = inputs['Area']
        qA = 0.5 * rho * Area * Vu ** 2
        outputs['Vd'] = Vd = Vu * (1 - 2 * a)
        outputs['Vr'] = 0.5 * (Vu + Vd)
        outputs['Ct'] = Ct = 4 * a * (1 - a)
        outputs['thrust'] = Ct * qA
        outputs['Cp'] = Cp = Ct * (1 - a)
        outputs['power'] = Cp * qA * Vu
        
    def compute_partials(self, inputs, J):
        a = inputs['a']
        Vu = inputs['Vu']
        Area = inputs['Area']
        rho = inputs['rho']

        a_times_area = a * Area
        one_minus_a = 1.0 - a
        a_area_rho_vu = a_times_area * rho * Vu
 
        J['Vr', 'a'] = -Vu
        J['Vr', 'Vu'] = one_minus_a
        J['Vd', 'a'] = -2.0 * Vu
        J['Ct', 'a'] = 4.0 - 8.0 * a
        J['thrust', 'a'] = 0.5 * rho * Vu**2 * Area * J['Ct', 'a']
        J['thrust', 'Area'] = 2.0 * Vu**2 * a * rho * one_minus_a
        J['thrust', 'Vu'] = 4.0 * a_area_rho_vu * one_minus_a
        J['Cp', 'a'] = 4.0 * a * (2.0 * a - 2.0) + 4.0 * one_minus_a**2
        J['power', 'a'] = 2.0 * Area * Vu**3 * a * rho * (2.0 * a - 2.0) + 2.0 * Area * Vu**3 * rho * one_minus_a**2
        J['power', 'Area'] = 2.0 * Vu**3 * a * rho * one_minus_a ** 2
        J['power', 'rho'] = 2.0 * a_times_area * Vu ** 3 * (one_minus_a)**2
        J['power', 'Vu'] = 6.0 * Area * Vu**2 * a * rho * one_minus_a**2<jupyter_output><empty_output><jupyter_text>The class declaration, `class ActuatorDisc(om.ExplicitComponent):` shows that our class, `ActuatorDisc` inherits off of the `ExplicitComponent` class in OpenMDAO.  In WISDEM, 99% of all coded components are of the `ExplicitComponent` class, so this is the most fundamental building block to get accustomed to.  Other types of components are described in the OpenMDAO docs [here](http://openmdao.org/twodocs/versions/latest/_srcdocs/packages/openmdao.components.html).The `ExplicitComponent` class provides a template for the user to: - Declare their input and output variables in the `setup` method- Calculate the outputs from the inputs in the `compute` method.  In an optimization loop, this is called at every iteration.- Calculate analytical gradients of outputs with respect to inputs in the `compute_partials` method.The variable declarations take the form of `self.add_input` or `self.add_output` where a variable name and default/initial vaue is assigned.  The value declaration also tells the OpenMDAO internals about the size and shape for any vector or multi-dimensional variables.  Other optional keywords that can help with code documentation and model consistency are `units=` and `desc=`. Working with analytical derivatives derivativesWe need to tell OpenMDAO which derivatives will need to be computed. That happens in the following lines:```pythonself.declare_partials('Vr', ['a', 'Vu'])self.declare_partials('Vd', 'a')self.declare_partials('Ct', 'a')self.declare_partials('thrust', ['a', 'Area', 'rho', 'Vu'])self.declare_partials('Cp', 'a')self.declare_partials('power', ['a', 'Area', 'rho', 'Vu'])```Note that lines like `self.declare_partials('Vr', ['a', 'Vu'])` references both the derivatives \\(\partial V_r / \partial a\\) and \\(\partial V_r / \partial V_u\\).The Jacobian in which we provide solutions to the derivatives is```pythonJ['Vr', 'a'] = -VuJ['Vr', 'Vu'] = one_minus_aJ['Vd', 'a'] = -2.0 * VuJ['Ct', 'a'] = 4.0 - 8.0 * aJ['thrust', 'a'] = 0.5 * rho * Vu**2 * Area * J['Ct', 'a']J['thrust', 'Area'] = 2.0 * Vu**2 * a * rho * one_minus_aJ['thrust', 'Vu'] = 4.0 * a_area_rho_vu * one_minus_aJ['Cp', 'a'] = 4.0 * a * (2.0 * a - 2.0) + 4.0 * one_minus_a**2J['power', 'a'] = 2.0 * Area * Vu**3 * a * rho * (2.0 * a - 2.0) + 2.0 * Area * Vu**3 * rho * one_minus_a**2J['power', 'Area'] = 2.0 * Vu**3 * a * rho * one_minus_a ** 2J['power', 'rho'] = 2.0 * a_times_area * Vu ** 3 * (one_minus_a)**2J['power', 'Vu'] = 6.0 * Area * Vu**2 * a * rho * one_minus_a**2``` In OpenMDAO, multiple components can be connected together inside of a Group.  There will be some other new elements to review, so let's take a look: Betz Group:<jupyter_code>class Betz(om.Group):
    """
    Group containing the actuator disc equations for deriving the Betz limit.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('a', 0.5)
        indeps.add_output('Area', 10.0, units='m**2')
        indeps.add_output('rho', 1.225, units='kg/m**3')
        indeps.add_output('Vu', 10.0, units='m/s')
        
        self.add_subsystem('a_disk', ActuatorDisc(), promotes=['a', 'Area', 'rho', 'Vu'])<jupyter_output><empty_output><jupyter_text>The `Betz` class derives off of the OpenMDAO `Group` class, which is typically the top-level class that is used in an analysis.  The OpenMDAO `Group` class allows you to cluster models in hierarchies. We can put multiple components in groups. We can also put other groups in groups.Components are added to groups with the `self.add_subsystem` command, which has two primary arguments.  The first is the string name to call the subsystem that is added and the second is the component or sub-group class instance.  A common optional argument is `promotes=`, which elevatest the input/ouput variable string names to the top-level namespace.  The `Betz` group shows examples where the `promotes=` can be passed a list of variable string names or the `'*'` wildcard to mean all input/output variables.The first subsystem that is added is an `IndepVarComp`, which are the independent variables of the problem.  Subsystem inputs that are not tied to other subsystem outputs should be connected to an independent variables.  For optimization problems, design variables must be part of an `IndepVarComp`.  In the Betz problem, we have `a`, `Area`, `rho`, and `Vu`.  Note that they are promoted to the top level namespace, otherwise we would have to access them by `'indeps.x'` and `'indeps.z'`.The next subsystem that is added is an instance of the component we created above:```pythonself.add_subsystem('a_disk', ActuatorDisc(), promotes=['a', 'Area', 'rho', 'Vu'])```The `promotes=` can also serve to connect variables.  In OpenMDAO, two variables with the same string name in the same namespace are automatically connected.  By promoting the same variable string names as in the `IndepCarComp`, they are automatically connected.  For variables that are not connected in this way, explicit connect statements are required, which is demonstrated in the next tutorial. Let's optimize our system!Even though we have all the pieces in a `Group`, we still need to put them into a `Problem` to be executed.  The `Problem` instance is where we can assign design variables, objective functions, and constraints.  It is also how the user interacts with the `Group` to set initial conditions and interrogate output values.First, we instantiate the `Problem` and assign an instance of `Betz` to be the root model:<jupyter_code>prob = om.Problem()
prob.model = Betz()<jupyter_output><empty_output><jupyter_text>Next we assign an optimization driver to the problem instance. If we only wanted to evaluate the model once and not optimize, then a driver is not needed:<jupyter_code>prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'<jupyter_output><empty_output><jupyter_text>With the optimization driver in place, we can assign design variables, objective(s), and constraints. Any IndepVarComp can be a design variable and any model output can be an objective or constraint.We want to maximize the objective, but OpenMDAO will want to minimize it as it is consistent with the standard optimization problem statement. So we minimize the negative to find the maximum. Note that `Cp` is not promoted from `a_disk`. Therefore we must reference it with `a_disk.Cp`<jupyter_code>prob.model.add_design_var('a', lower=0.0, upper=1.0)
prob.model.add_design_var('Area', lower=0.0, upper=1.0)
prob.model.add_objective('a_disk.Cp', scaler=-1.0)<jupyter_output><empty_output><jupyter_text>Now we can run the optimization:<jupyter_code>prob.setup()
prob.run_driver()<jupyter_output>Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.5925925906659251
            Iterations: 5
            Function evaluations: 6
            Gradient evaluations: 5
Optimization Complete
-----------------------------------<jupyter_text>Finally, the result:Above, we see a summary of the steps in our optimization. Don't worry about the output `False` for now. Next, we print out the values we care about and list all of the inputs and outputs that are problem used.<jupyter_code>print('Coefficient of power Cp = ', prob['a_disk.Cp'])
print('Induction factor a =', prob['a'])
print('Rotor disc Area =', prob['Area'], 'm^2')
all_inputs  = prob.model.list_inputs(values=True)
all_outputs = prob.model.list_outputs(values=True)<jupyter_output>Coefficient of power Cp =  [0.59259259]
Induction factor a = [0.33335528]
Rotor disc Area = [1.] m^2
4 Input(s) in 'model'
---------------------

varname   value       
--------  ------------
top
  a_disk
    a     [0.33335528]
    Area  [1.]        
    rho   [1.225]     
    Vu    [10.]       


10 Explicit Output(s) in 'model'
--------------------------------

varname     value         
----------  --------------
top
  indeps
    a       [0.33335528]  
    Area    [1.]          
    rho     [1.225]       
    Vu      [10.]         
  a_disk
    Vr      [6.6664472]   
    Vd      [3.33289439]  
    Ct      [0.88891815]  
    Cp      [0.59259259]  
    power   [362.96296178]
    thrust  [54.44623668] 


0 Implicit Output(s) in 'model'
-------------------------------
{'fraction_non_alphanumeric': 0.1331332796603967, 'fraction_numerical': 0.022323062285003293, 'mean_word_length': 3.520013231888852, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '156470', 'license': 'Apache-2.0', 'n_tokens_mistral': 4537, 'n_tokens_neox': 4220, 'n_words': 1579, 'path': 'tutorial-notebooks/02_betz_limit.ipynb', 'repo_name': 'dzalkind/WISDEM', 'size': '13663'}