Cantilevered Beam

Note

This problem can be resolved without the surrogate modeling functionality. Indeed, this functionality is useful when the objective function is computationally expensive. The purpose of this example is to illustrate the use of an array external function and the use of surrogate modeling on a simple and computationally inexpensive problem.

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

../_images/cantilevered_beam.svg

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, an O_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.

Execution:
localsolver cantilevered_beam.lsp [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=%LS_HOME%\bin\python
python cantilevered_beam.py
Execution (Linux)
export PYTHONPATH=/opt/localsolver_12_0/bin/python
python cantilevered_beam.py
import localsolver
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 localsolver.LocalSolver() as ls:
        # Declare the optimization model
        model = ls.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 solver
        surrogate_params.evaluation_limit = evaluation_limit

        ls.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%LS_HOME%\include /link %LS_HOME%\bin\localsolver120.lib
cantilevered_beam
Compilation / Execution (Linux)
g++ cantilevered_beam.cpp -I/opt/localsolver_12_0/include -llocalsolver120 -lpthread -o cantilevered_beam
cantilevered_beam
#include "localsolver.h"
#include <vector>
#include <iostream>
#include <fstream>

using namespace std;
using namespace localsolver;

/* 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 LSArrayExternalFunction<lsdouble> {

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

        // Big computation
        lsdouble 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
        lsdouble g1 = P * L * fH / (2 * I);
        lsdouble g2 = P * pow(L, 3) / (3 * E * I);

        // Objective function computations
        lsdouble 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:
    // LocalSolver
    LocalSolver localsolver;

    // LS Program variables
    LSExpression H;
    LSExpression h1;
    LSExpression b1;
    LSExpression b2;

    LSExpression objective;

    void solve(int evaluationLimit){
        // Declare the optimization model
        LSModel model = localsolver.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;
        LSExpression func = model.createExternalFunction(&evaluate);
        LSExpression results = func(H, h1, b1, b2);

        // Enable surrogate modeling
        LSExternalContext context = func.getExternalContext();
        LSSurrogateParameters 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 solver
        surrogateParams.setEvaluationLimit(evaluationLimit);

        localsolver.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 %LS_HOME%\bin\localsolvernet.dll .
csc CantileveredBeam.cs /reference:localsolvernet.dll
CantileveredBeam
using System;
using System.IO;
using localsolver;

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(LSExternalArgumentValues 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};
        }
    }

    // LocalSolver
    private LocalSolver localsolver;

    // LS Program variables
    private LSExpression H;
    private LSExpression h1;
    private LSExpression b1;
    private LSExpression b2;

    private LSExpression objective;

    public CantileveredBeam() {
        localsolver = new LocalSolver();
    }

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

    public void Solve(int evaluationLimit) {
        // Declare the optimization model
        LSModel model = localsolver.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();
        LSDoubleArrayExternalFunction func = new LSDoubleArrayExternalFunction(evaluateFunction.Call);
        LSExpression funcExpr = model.DoubleArrayExternalFunction(func);
        LSExpression results = model.Call(funcExpr, H, h1, b1, b2);

        // Enable surrogate modeling
        LSExternalContext context = funcExpr.GetExternalContext();
        LSSurrogateParameters 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 solver
        surrogateParams.SetEvaluationLimit(evaluationLimit);
        localsolver.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 %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. CantileveredBeam
Compilation / Execution (Linux)
javac CantileveredBeam.java -cp /opt/localsolver_12_0/bin/localsolver.jar
java -cp /opt/localsolver_12_0/bin/localsolver.jar:. CantileveredBeam
import java.io.*;
import localsolver.*;

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 LSDoubleArrayExternalFunction {
        @Override
        public double[] call(LSExternalArgumentValues 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};
        }
    }

    // LocalSolver
    private final LocalSolver localsolver;

    // LS Program variables
    private LSExpression H;
    private LSExpression h1;
    private LSExpression b1;
    private LSExpression b2;

    private LSExpression objective;

    private CantileveredBeam(LocalSolver localsolver) {
        this.localsolver = localsolver;
    }

    private void solve(int evaluationLimit) {
        // Declare the optimization model
        LSModel model = localsolver.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
        LSExpression func = model.doubleArrayExternalFunction(new Evaluate());
        LSExpression results = model.call(func, H, h1, b1, b2);

        // Enable surrogate modeling
        LSSurrogateParameters 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 solver
        surrogateParams.setEvaluationLimit(evaluationLimit);

        localsolver.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 (LocalSolver localsolver = new LocalSolver()) {
            CantileveredBeam model = new CantileveredBeam(localsolver);
            model.solve(Integer.parseInt(strEvaluationLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}