Optimal Bucket Problem

Problem

What is the optimal shape for a bucket? In the Optimal Bucket Problem, we wish to design a bucket that maximizes the volume of fluid it can hold while respecting the available surface material. A bucket is defined by three values: the radius of its bottom disc, the radius of its top opening, and its height, respectively noted r, R, and h. The problem consists in choosing the values of r, R, and h so as to maximize the bucket’s volume under the constraint that its surface does not exceed that of the available material.

For more details, see DataGenetics.

Principles learned

Program

The available material to build the bucket is a flat disk with a radius equal to 1. Its surface area is S=π. Without using more material than this amount, we try to build a bucket that holds the largest volume.

The model has three float decision variables. They represent the three quantities defining the bucket’s shape: the radius of its bottom disk r, the radius of its top opening R, and its height h. The surface and the volume of the bucket are completely determined from these three quantities. Therefore, they do not need to be decision variables but intermediate expressions. The surface of the bucket is given by the expression S = π*r² + π(R+r)sqrt((R-r)²+h²). After computing it, we constrain it to be lower than π. We can then compute the volume of the bucket, given by the expression V = (π*h)/3 * (R²+Rr+r²), and we maximize it.

Execution
hexaly optimal_bucket.hxm [hxTimeLimit=] [solFileName=]
use io;

/* Declare the optimization model */
function model() {
    PI = 3.14159265359;

    // Numerical decisions
    R <- float(0, 1);
    r <- float(0, 1);
    h <- float(0, 1);

    // Surface must not exceed the surface of the plain disc
    surface <- PI * pow(r, 2) + PI * (R + r) * sqrt(pow(R - r, 2) + pow(h, 2));
    constraint surface <= PI;

    // Maximize the volume
    volume <- PI * h / 3 * (pow(R, 2) + R * r + pow(r, 2));
    maximize volume;
}

/* Parametrize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 2;
}

/* Write the solution in a file with the following format:
 *  - surface and volume of the bucket
 *  - values of R, r and h */
function output() {
    if (solFileName == nil) return;
    local solFile = io.openWrite(solFileName);
    solFile.println(surface.value, "  ", volume.value);
    solFile.println(R.value, "  ", r.value, "  ", h.value);
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python optimal_bucket.py
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python optimal_bucket.py 
import hexaly.optimizer
import sys

with hexaly.optimizer.HexalyOptimizer() as optimizer:
    PI = 3.14159265359

    #
    # Declare the optimization model
    #
    m = optimizer.model

    # Numerical decisions
    R = m.float(0, 1)
    r = m.float(0, 1)
    h = m.float(0, 1)

    # Surface must not exceed the surface of the plain disc
    surface = PI * r ** 2 + PI * (R + r) * m.sqrt((R - r) ** 2 + h ** 2)
    m.constraint(surface <= PI)

    # Maximize the volume
    volume = PI * h / 3 * (R ** 2 + R * r + r ** 2)
    m.maximize(volume)

    m.close()

    #
    # Parametrize the optimizer
    #
    if len(sys.argv) >= 3:
        optimizer.param.time_limit = int(sys.argv[2])
    else:
        optimizer.param.time_limit = 2

    optimizer.solve()

    #
    # Write the solution in a file with the following format:
    #  - surface and volume of the bucket
    #  - values of R, r and h
    #
    if len(sys.argv) >= 2:
        with open(sys.argv[1], 'w') as f:
            f.write("%f %f\n" % (surface.value, volume.value))
            f.write("%f %f %f\n" % (R.value, r.value, h.value))
Compilation / Execution (Windows)
cl /EHsc optimal_bucket.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
optimal_bucket
Compilation / Execution (Linux)
g++ optimal_bucket.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o optimal_bucket
./optimal_bucket
#include "optimizer/hexalyoptimizer.h"
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

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

    // Hexaly Program variables
    HxExpression R;
    HxExpression r;
    HxExpression h;

    HxExpression surface;
    HxExpression volume;

    void solve(int limit) {
        hxdouble PI = 3.14159265359;

        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Numerical decisions
        R = model.floatVar(0.0, 1.0);
        r = model.floatVar(0.0, 1.0);
        h = model.floatVar(0.0, 1.0);

        // Surface must not exceed the surface of the plain disc
        surface = PI * model.pow(r, 2) + PI * (R + r) * model.sqrt(model.pow(R - r, 2) + model.pow(h, 2));
        model.constraint(model.leq(surface, PI));

        // Maximize the volume
        volume = PI * h / 3 * (model.pow(R, 2) + R * r + model.pow(r, 2));
        model.maximize(volume);

        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     *  - surface and volume of the bucket
     *  - values of R, r and h */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << surface.getDoubleValue() << " " << volume.getDoubleValue() << endl;
        outfile << R.getDoubleValue() << " " << r.getDoubleValue() << " " << h.getDoubleValue() << endl;
    }
};

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

    try {
        OptimalBucket model;
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred:" << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc OptimalBucket.cs /reference:Hexaly.NET.dll
OptimalBucket
using System;
using System.IO;
using Hexaly.Optimizer;

public class OptimalBucket : IDisposable
{
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Hexaly Program variables
    HxExpression R;
    HxExpression r;
    HxExpression h;

    HxExpression surface;
    HxExpression volume;

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

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

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

        // Numerical decisions
        R = model.Float(0, 1);
        r = model.Float(0, 1);
        h = model.Float(0, 1);

        // Surface must not exceed the surface of the plain disc
        surface =
            Math.PI * model.Pow(r, 2)
            + Math.PI * (R + r) * model.Sqrt(model.Pow(R - r, 2) + model.Pow(h, 2));
        model.AddConstraint(surface <= Math.PI);

        // Maximize the volume
        volume = Math.PI * h / 3 * (model.Pow(R, 2) + R * r + model.Pow(r, 2));
        model.Maximize(volume);

        model.Close();

        // Parametrize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

    // Write the solution in a file with the following format:
    //  - surface and volume of the bucket
    //  - values of R, r and h
    public void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(surface.GetDoubleValue() + " " + volume.GetDoubleValue());
            output.WriteLine(
                R.GetDoubleValue() + " " + r.GetDoubleValue() + " " + h.GetDoubleValue()
            );
        }
    }

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

        using (OptimalBucket model = new OptimalBucket())
        {
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac OptimalBucket.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. OptimalBucket
Compilation / Execution (Linux)
javac OptimalBucket.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. OptimalBucket
import java.io.*;
import com.hexaly.optimizer.*;

public class OptimalBucket {
    private static final double PI = 3.14159265359;

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Hexaly Program variables
    private HxExpression R;
    private HxExpression r;
    private HxExpression h;

    private HxExpression surface;
    private HxExpression volume;

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

    private void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();
        HxExpression piConst = model.createConstant(PI);

        // Numerical decisions
        R = model.floatVar(0, 1);
        r = model.floatVar(0, 1);
        h = model.floatVar(0, 1);

        // surface = PI*r^2 + PI*(R+r)*sqrt((R - r)^2 + h^2)
        HxExpression s1 = model.prod(piConst, r, r);
        HxExpression s2 = model.pow(model.sub(R, r), 2);
        HxExpression s3 = model.pow(h, 2);
        HxExpression s4 = model.sqrt(model.sum(s2, s3));
        HxExpression s5 = model.sum(R, r);
        HxExpression s6 = model.prod(piConst, s5, s4);
        surface = model.sum(s1, s6);

        // Surface must not exceed the surface of the plain disc
        model.addConstraint(model.leq(surface, PI));

        HxExpression v1 = model.pow(R, 2);
        HxExpression v2 = model.prod(R, r);
        HxExpression v3 = model.pow(r, 2);

        // volume = PI*h/3*(R^2 + R*r + r^2)
        volume = model.prod(piConst, model.div(h, 3), model.sum(v1, v2, v3));

        // Maximize the volume
        model.maximize(volume);

        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     * - surface and volume of the bucket
     * - values of R, r and h */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(surface.getDoubleValue() + " " + volume.getDoubleValue());
            output.println(R.getDoubleValue() + " " + r.getDoubleValue() + " " + h.getDoubleValue());
        }
    }

    public static void main(String[] args) {

        String outputFile = args.length > 0 ? args[0] : null;
        String strTimeLimit = args.length > 1 ? args[1] : "2";

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

Discover the ease of use and performance of Hexaly