Cantilevered Beam Problem

Problem

The Cantilevered Beam Problem consists in designing the cross-sectional of the I-beam in order to get the minimum volume without deforming or breaking the beam under a certain stress. For more details, see Wikipedia. We note L, H, h1, b1, and b2 the dimensions of the beam. They are defined according to the following picture.

The volume of the beam is given by the equation:

V = f (H, h_{1}, b_{1}, b_{2}) = [2 * h_{1} * b_{1} + (H - 2 * h_{1}) * b_{2}] * L

The constraints are:

The maximum bending stress at the root of the beam, defined by:

g 1 (H, h_{1}, b_{1}, b_{2}) = (P * L * H) / (2 * I)

The maximum deflection at the tip of the beam, defined by:

g 2 (H, h_{1}, b_{1}, b_{2}) = (P * L^{3}) / (3 * E * I)

Principles learned

Create an external function returning multiple values
Enable surrogate modeling on an external function
Set an evaluation limit to the function

Note: The purpose of this example is to illustrate how to define an external function returning an array, and how to use surrogate modeling on a simple problem. Since the external function used in this problem is computationally inexpensive, it can be solved without the surrogate modeling functionality. Indeed, surrogate modeling is only useful when the objective function is computationally expensive. In the example, we use it for illustration purposes only.

Program

The Hexaly model for the Cantilevered Beam Problem uses int and float decision variables. We start by declaring the three float decision variables H, b1, and b2. The domains of these variables are respectively [3.0, 7.0], [2.0, 12.0], and [3.0, 7.0]. The last decision variable h1 is discrete. It is declared as an integer representing the index in the array containing the possible values, which are {0.1, 0.26, 0.35, 0.5, 0.65, 0.75, 0.9, 1.0}. The domain of this variable is then [0, 7].

The volume, bending stress, and deflection of the beam are computed using an external function. It receives the argument values via an HxExternalArgumentValues object, and returns an array containing the values of the three functions for these beam characteristics. To compute the values returned by the external function, a call expression is created. We use the first two return values (bending stress and deflection) to constrain the model. The third return value (volume) is the objective function to be minimized.

To use the surrogate modeling feature, we call the enableSurrogateModeling method available on the HxExternalContext of the function. This method returns the HxSurrogateParameters, on which the maximum number of calls to the function can be set. As the function is usually computationally expensive in practice, it is very useful to limit the search to a reasonable time.

Hexaly Modeler
Python
C++
C#
Java

Download

Execution: hexaly cantilevered_beam.hxm [evaluationLimit=] [solFileName=]

use io;

/* Constant declaration */
function input() {
    P = 1000;
    E = 10000000;
    L = 36;
    possibleValues = {0.1, 0.25, 0.35, 0.5, 0.65, 0.75, 0.9, 1.0};
}

/* External function */
function evaluate(args){
    fH = args[0];
    fh1 = possibleValues[args[1]];
    fb1 = args[2];
    fb2 = args[3];

    // Big computation
    I = 1.0 / 12.0 * fb2 * pow(fH - 2 * fh1, 3) + 2 * (1.0 / 12.0 * fb1
        * pow(fh1, 3) + fb1 * fh1 * pow(fH - fh1, 2) / 4.0);

    // Constraint computations
    g1 = P * L * fH / (2 * I);
    g2 = P * pow(L, 3) / (3 * E * I);

    // Objective function computation
    f = (2 * fh1 * fb1 + (fH - 2 * fh1) * fb2) * L;

    return {g1, g2, f};
}

/* Declare the optimization model */
function model() {
    // Numerical decisions
    H <- float(3.0, 7.0);
    h1 <- int(0, 7);
    b1 <- float(2.0, 12.0);
    b2 <- float(0.1, 2.0);

    // Create and call the external function
    func <- doubleArrayExternalFunction(evaluate);
    results <- call(func, H, h1, b1, b2);

    // Enable surrogate modeling
    surrogateParams = func.context.enableSurrogateModeling();

    // Constraint on bending stress
    constraint results[0] <= 5000;
    // Constraint on deflection at the tip of the beam
    constraint results[1] <= 0.10;

    objective <- results[2];
    minimize objective;
}

/* Parameterize the solver */
function param() {
    if (evaluationLimit == nil) surrogateParams.evaluationLimit = 30;
    else surrogateParams.evaluationLimit = evaluationLimit;
}

/* Write the solution in a file with the following format:
 * - The value of the minimum found
 * - The location (H; h1; b1; b2) of the minimum
*/
function output() {
    if (solFileName != nil) {
        local solFile = io.openWrite(solFileName);
        solFile.println("Objective value: ", objective.value);
        solFile.println("Point (H;h1;b1;b2): (" + H.value + ";"
            + possibleValues[h1.value] + ";" + b1.value + ";" + b2.value + ")");
    }
}

Execution (Windows): set PYTHONPATH=%HX_HOME%\bin\python

python cantilevered_beam.py

Execution (Linux): export PYTHONPATH=/opt/hexaly_13_0/bin/python

python cantilevered_beam.py

import hexaly.optimizer
import sys

# Constant declaration
P = 1000
E = 10.0e6
L = 36
possibleValues = [0.1, 0.26, 0.35, 0.5, 0.65, 0.75, 0.9, 1.0]

# External function
def evaluate(arguments_values):
    # Argument retrieval
    fH = arguments_values[0]
    fh1 = possibleValues[arguments_values[1]]
    fb1 = arguments_values[2]
    fb2 = arguments_values[3]

    # Big computation
    I = 1.0 / 12.0 * fb2 * pow(fH - 2 * fh1, 3) + 2 * (1.0 / 12.0 * fb1
        * pow(fh1, 3) + fb1 * fh1 * pow(fH - fh1, 2) / 4.0)

    # Constraint computations
    g1 = P * L * fH / (2 * I)
    g2 = P * L**3 / (3 * E * I)

    # Objective function computation
    f = (2 * fh1 * fb1 + (fH - 2 * fh1) * fb2) * L

    return g1, g2, f

def main(evaluation_limit, output_file):
    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        # Declare the optimization model
        model = optimizer.model

        # Numerical decisions
        H = model.float(3.0, 7.0)
        h1 = model.int(0, 7)
        b1 = model.float(2.0, 12.0)
        b2 = model.float(0.1, 2.0)

        # Create and call the external function
        func = model.create_double_array_external_function(evaluate)
        func_call = model.call(func)
        # Add the operands
        func_call.add_operand(H)
        func_call.add_operand(h1)
        func_call.add_operand(b1)
        func_call.add_operand(b2)

        # Enable surrogate modeling
        surrogate_params = func.external_context.enable_surrogate_modeling()

        # Constraint on bending stress
        model.constraint(func_call[0] <= 5000)
        # Constraint on deflection at the tip of the beam
        model.constraint(func_call[1] <= 0.10)

        objective = func_call[2]
        model.minimize(objective)
        model.close()

        # Parameterize the optimizer
        surrogate_params.evaluation_limit = evaluation_limit

        optimizer.solve()

        # Write the solution in a file with the following format:
        # - The value of the minimum found
        # - The location (H; h1; b1; b2) of the minimum
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("Objective value: " + str(objective.value) + "\n")
                f.write("Point (H;h1;b1;b2): (" + str(H.value) + ";"
                    + str(possibleValues[h1.value]) + ";" + str(b1.value) + ";"
                    + str(b2.value) + ")")


if __name__ == '__main__':
    output_file = sys.argv[1] if len(sys.argv) > 1 else None
    evaluation_limit = int(sys.argv[2]) if len(sys.argv) > 2 else 30

    main(evaluation_limit, output_file)

Compilation / Execution (Windows): cl /EHsc cantilevered_beam.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib

cantilevered_beam

Compilation / Execution (Linux): g++ cantilevered_beam.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o cantilevered_beam

./cantilevered_beam

#include "optimizer/hexalyoptimizer.h"
#include <vector>
#include <iostream>
#include <fstream>

using namespace std;
using namespace hexaly;

/* Constant declaration */
const int P = 1000;
const int E = 10000000;
const int L = 36;
const float possibleValues[8] = {0.1, 0.25, 0.35, 0.5, 0.65, 0.75, 0.9, 1.0};

/* External function */
class Evaluate : public HxArrayExternalFunction<hxdouble> {

    void call(const HxExternalArgumentValues& argumentValues,
                vector<hxdouble>& resultValues) override {
        // Argument retrieval
        hxdouble fH = argumentValues.getDoubleValue(0);
        hxdouble fh1 = possibleValues[argumentValues.getIntValue(1)];
        hxdouble fb1 = argumentValues.getDoubleValue(2);
        hxdouble fb2 = argumentValues.getDoubleValue(3);

        // Big computation
        hxdouble I = 1.0 / 12.0 * fb2 * pow(fH - 2 * fh1, 3) + 2 * (1.0 / 12.0
            * fb1 * pow(fh1, 3) + fb1 * fh1 * pow(fH - fh1, 2) / 4.0);

        // Constraint computations
        hxdouble g1 = P * L * fH / (2 * I);
        hxdouble g2 = P * pow(L, 3) / (3 * E * I);

        // Objective function computations
        hxdouble f = (2 * fh1 * fb1 + (fH - 2 * fh1) * fb2) * L;

        // Return results
        resultValues.clear();
        resultValues.push_back(g1);
        resultValues.push_back(g2);
        resultValues.push_back(f);
    }
};

class CantileveredBeam {
public:
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Hexaly Program variables
    HxExpression H;
    HxExpression h1;
    HxExpression b1;
    HxExpression b2;

    HxExpression objective;

    void solve(int evaluationLimit){
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Numerical decisions
        H = model.floatVar(3.0, 7.0);
        h1 = model.intVar(0, 7);
        b1 = model.floatVar(2.0, 12.0);
        b2 = model.floatVar(0.1, 2.0);

        // Create and call the external function
        Evaluate evaluate;
        HxExpression func = model.createExternalFunction(&evaluate);
        HxExpression results = func(H, h1, b1, b2);

        // Enable surrogate modeling
        HxExternalContext context = func.getExternalContext();
        HxSurrogateParameters surrogateParams = context.enableSurrogateModeling();

        // Constraint on bending stress
        model.constraint(results[0] <= 5000);
        // Constraint on deflection at the tip of the beam
        model.constraint(results[1] <= 0.10);

        objective = results[2];
        model.minimize(objective);
        model.close();

        // Parameterize the optimizer
        surrogateParams.setEvaluationLimit(evaluationLimit);

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
    * - The value of the minimum found
    * - The location (H; h1; b1; b2) of the minimum
    */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());
        outfile << "Objective value: " << objective.getDoubleValue() << endl;
        outfile << "Point (H;h1;b1;b2): (" << H.getDoubleValue() << ";"
            << possibleValues[h1.getIntValue()] << ";" << b1.getDoubleValue()
            << ";" << b2.getDoubleValue() << ")";
    }
};

int main(int argc, char** argv){
    const char* solFile = argc > 1 ? argv[1] : NULL;
    const char* strEvaluationLimit = argc > 2 ? argv[2] : "30";

    try {
        CantileveredBeam model;
        model.solve(atoi(strEvaluationLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
    return 0;
}

Compilation / Execution (Windows): copy %HX_HOME%\bin\Hexaly.NET.dll .

csc CantileveredBeam.cs /reference:Hexaly.NET.dll

CantileveredBeam

using System;
using System.IO;
using Hexaly.Optimizer;

public class CantileveredBeam : IDisposable
{
    public class Evaluate {
        /* Constant declaration */
        const int P = 1000;
        const int E = 10000000;
        const int L = 36;
        public static double[] possibleValues = {0.1, 0.25, 0.35, 0.5, 0.65, 0.75, 0.9, 1.0};

        /* External function */
        public double[] Call(HxExternalArgumentValues argumentValues) {
            // Argument retrieval
            double fH = argumentValues.GetDoubleValue(0);
            double fh1 = possibleValues[(int)argumentValues.GetIntValue(1)];
            double fb1 = argumentValues.GetDoubleValue(2);
            double fb2 = argumentValues.GetDoubleValue(3);

            // Big computation
            double I = 1.0 / 12.0 * fb2 * Math.Pow(fH - 2 * fh1, 3) + 2 * (1.0 / 12.0
                * fb1 * Math.Pow(fh1, 3) + fb1 * fh1 * Math.Pow(fH - fh1, 2) / 4.0);

            // Constraint computations
            double g1 = P * L * fH / (2 * I);
            double g2 = P * Math.Pow(L, 3) / (3* E * I);

            // Objective function computation
            double f = (2 * fh1 * fb1 + (fH - 2 * fh1) * fb2) * L;

            return new double[]{g1, g2, f};
        }
    }

    // Hexaly Optimizer
    private HexalyOptimizer optimizer;

    // Hexaly Program variables
    private HxExpression H;
    private HxExpression h1;
    private HxExpression b1;
    private HxExpression b2;

    private HxExpression objective;

    public CantileveredBeam() {
        optimizer = new HexalyOptimizer();
    }

    public void Dispose() {
        if (optimizer != null)
            optimizer.Dispose();
    }

    public void Solve(int evaluationLimit) {
        // Declare the optimization model
        HxModel model = optimizer.GetModel();

        // Numerical decisions
        H = model.Float(3.0, 7.0);
        h1 = model.Int(0, 7);
        b1 = model.Float(2.0, 12.0);
        b2 = model.Float(0.1, 2.0);

        // Create and call the external function
        Evaluate evaluateFunction = new Evaluate();
        HxDoubleArrayExternalFunction func = new HxDoubleArrayExternalFunction(evaluateFunction.Call);
        HxExpression funcExpr = model.DoubleArrayExternalFunction(func);
        HxExpression results = model.Call(funcExpr, H, h1, b1, b2);

        // Enable surrogate modeling
        HxExternalContext context = funcExpr.GetExternalContext();
        HxSurrogateParameters surrogateParams = context.EnableSurrogateModeling();

        // Constraint on bending stress
        model.Constraint(results[0] <= 5000);
        // Constraint on deflection at the tip of the beam
        model.Constraint(results[1] <= 0.10);

        objective = results[2];
        model.Minimize(objective);
        model.Close();

        // Parameterize the optimizer
        surrogateParams.SetEvaluationLimit(evaluationLimit);
        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
    * - The value of the minimum found
    * - The location (H; h1; b1; b2) of the minimum
    */
    public void WriteSolution(string fileName) {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("Objective value: " + objective.GetDoubleValue());
            output.WriteLine("Point (H;h1;b1;b2): (" + H.GetDoubleValue()
                + ";" + Evaluate.possibleValues[(int)h1.GetIntValue()] + ";" + b1.GetDoubleValue()
                + ";" + b2.GetDoubleValue() + ")");
        }
    }

    public static void Main(string[] args)
    {
        string outputFile = args.Length > 0 ? args[0] : null;
        string strEvaluationLimit = args.Length > 1 ? args[1] : "30";

        using (CantileveredBeam model = new CantileveredBeam())
        {
            model.Solve(int.Parse(strEvaluationLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}

Compilation / Execution (Windows): javac CantileveredBeam.java -cp %HX_HOME%\bin\hexaly.jar

java -cp %HX_HOME%\bin\hexaly.jar;. CantileveredBeam

Compilation / Execution (Linux): javac CantileveredBeam.java -cp /opt/hexaly_13_0/bin/hexaly.jar

java -cp /opt/hexaly_13_0/bin/hexaly.jar:. CantileveredBeam

import java.io.*;
import com.hexaly.optimizer.*;

public class CantileveredBeam {

    /* Constant declaration */
    final static int P = 1000;
    final static int E = 10000000;
    final static int L = 36;
    final static double[] possibleValues = {0.1, 0.25, 0.35, 0.5, 0.65, 0.75, 0.9, 1.0};

    /* External function */
    private static class Evaluate implements HxDoubleArrayExternalFunction {
        @Override
        public double[] call(HxExternalArgumentValues argumentValues){
            // Argument retrieval
            double fH = argumentValues.getDoubleValue(0);
            double fh1 = possibleValues[(int)argumentValues.getIntValue(1)];
            double fb1 = argumentValues.getDoubleValue(2);
            double fb2 = argumentValues.getDoubleValue(3);

            // Big computation
            double I = 1.0 / 12.0 * fb2 * Math.pow(fH - 2 * fh1, 3) + 2 * (1.0 / 12.0
                * fb1 * Math.pow(fh1, 3) + fb1 * fh1 * Math.pow(fH - fh1, 2) / 4.0);

            // Constraint computations
            double g1 = P * L * fH / (2 * I);
            double g2 = P * Math.pow(L, 3) / (3* E * I);

            // Objective function computation
            double f = (2 * fh1 * fb1 + (fH - 2 * fh1) * fb2) * L;

            return new double[]{g1, g2, f};
        }
    }

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Hexaly Program variables
    private HxExpression H;
    private HxExpression h1;
    private HxExpression b1;
    private HxExpression b2;

    private HxExpression objective;

    private CantileveredBeam(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

    private void solve(int evaluationLimit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Numerical decisions
        H = model.floatVar(3.0, 7.0);
        h1 = model.intVar(0, 7);
        b1 = model.floatVar(2.0, 12.0);
        b2 = model.floatVar(0.1, 2.0);

        // Create and call the external function
        HxExpression func = model.doubleArrayExternalFunction(new Evaluate());
        HxExpression results = model.call(func, H, h1, b1, b2);

        // Enable surrogate modeling
        HxSurrogateParameters surrogateParams = func.getExternalContext().enableSurrogateModeling();

        // Constraint on bending stress
        model.constraint(model.leq(model.at(results, 0), 5000));
        // Constraint on deflection at the tip of the beam
        model.constraint(model.leq(model.at(results, 1), 0.10));

        objective = model.at(results, 2);
        model.minimize(objective);
        model.close();

        // Parameterize the optimizer
        surrogateParams.setEvaluationLimit(evaluationLimit);

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
    * - The value of the minimum found
    * - The location (H; h1; b1; b2) of the minimum
    */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println("Objective value: " + objective.getDoubleValue());
            output.println("Point (H;h1;b1;b2): (" + H.getDoubleValue()
                + ";" + possibleValues[(int)h1.getIntValue()] + ";" + b1.getDoubleValue()
                + ";" + b2.getDoubleValue() + ")");
        }
    }

    public static void main(String[] args) {
        String outputFile = args.length > 0 ? args[0] : null;
        String strEvaluationLimit = args.length > 1 ? args[1] : "30";

        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            CantileveredBeam model = new CantileveredBeam(optimizer);
            model.solve(Integer.parseInt(strEvaluationLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}

Problem

Principles learned

Program

Similar examples