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)\]
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.
- 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\pythonpython cantilevered_beam.py
- Execution (Linux)
- export PYTHONPATH=/opt/localsolver_12_5/bin/pythonpython 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\localsolver125.libcantilevered_beam
- Compilation / Execution (Linux)
- g++ cantilevered_beam.cpp -I/opt/localsolver_12_5/include -llocalsolver125 -lpthread -o cantilevered_beamcantilevered_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.dllCantileveredBeam
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.jarjava -cp %LS_HOME%\bin\localsolver.jar;. CantileveredBeam
- Compilation / Execution (Linux)
- javac CantileveredBeam.java -cp /opt/localsolver_12_5/bin/localsolver.jarjava -cp /opt/localsolver_12_5/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);
}
}
}