Cantilevered Beam¶

Note

The purpose of this example is to illustrate the use of an array external function and the use of surrogate modeling on a simple problem.

Since, in this problem, the external function is computationally inexpensive, it can be solved without the surrogate modeling functionality. Indeed, surrogate modeling is useful when the objective function is computationally expensive. In the example, it is used for illustration purpose only.

Principles learned¶

Create an external function returning multiple values
Enable the surrogate modeling on an external function

Set an evaluation limit to the function

Problem¶

How to minimize the volume of a cantilevered I-beam subjected to a tip load?

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.

The volume of the beam is defined by:

\[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:

\[g1(H, h_1, b_1, b_2) = (P * L * H) / (2 * I)\]
The maximum deflection at the tip of the beam:

\[g2(H, h_1, b_1, b_2) = (P * L^3) / (3 * E * I)\]

Download the example

Program¶

The objective and the constraint functions are defined by an external function. It receives the argument values via the LSExternalArgumentValues, and returns the value of the functions at this point.

Three floating decision variables H, b1 and b2 are declared. The domains of these variables are respectively [3.0, 7.0], [2.0, 12.0] and [3.0, 7.0]. The final decision variable, h1 is discrete. It is declared as an integer representing the index of 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 index is [0, 7].

To calculate the values returned by the external function, a call expression is created. The two first returned expressions will be used to constraint the model. The third one will be minimized.

To use the surrogate modeling feature, the method enableSurrogateModeling available on the LSExternalContext of the function is called. This method returns the LSSurrogateParameters, on which the maximum number of calls to the function can be set, as the function is usually computationally expensive.

Hexaly Modeler
Python
C++
C#
Java

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(fH, indexH1, fb1, fb2){
    fh1 = possibleValues[indexH1];

    // 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);
        }
    }
}