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:
The constraints are:
- The maximum bending stress at the root of the beam, defined by:
- The maximum deflection at the tip of the beam, defined by:
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.
- 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\pythonpython cantilevered_beam.py
- Execution (Linux)
-
export PYTHONPATH=/opt/hexaly_13_0/bin/pythonpython 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.libcantilevered_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.dllCantileveredBeam
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.jarjava -cp %HX_HOME%\bin\hexaly.jar;. CantileveredBeam
- Compilation / Execution (Linux)
-
javac CantileveredBeam.java -cp /opt/hexaly_13_0/bin/hexaly.jarjava -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);
}
}
}