借助 chaptgpt 和 deepseek，成功实现了c++上的多阶段报童模型的动态规划。花费了几天，将以前的 java 程序用 c++ 实现。

总结：

• c++ 还是比 java 快点，30个阶段快了零点几秒
• c++ 使用了 unordered_map ，存储递归数据
• java 使用了 ConcurrentSkipListMap 存储递归数据，这个可以按照排序器自动排序
• 若 c++ 也用可以排序的 map，速度反而比 java 慢了。理论上c++会快，但估计需要其他的一些功能设置
• c++ 运行时要开启 -o2 或 -o3 优化加速

C++ 代码

//
// Created by Zhen Chen on 2025/2/26.
//

#include <chrono>
#include <iostream>
#include <limits>
#include <boost/functional/hash.hpp>
#include <unordered_map>
#include <iostream>
#include <map>
#include <span>
#include <csignal>

class State {
    int period{}; // c++11, {} 值初始化，默认为 0
    double initialInventory{};

public:
    State();

    explicit State(int period, double initialInventory);

    [[nodiscard]] double getInitialInventory() const;

    [[nodiscard]] int getPeriod() const;

    void print() const;

    // hashmap must define operator == and a struct to compute hash
    bool operator==(const State &other) const {
        // 需要定义 `==`
        // const MyClass &other	保证 other 参数不可修改
        // const 在函数结尾 保证当前对象(this) 不可修改
        // 不会修改成员变量的方法 都可以在函数声明的结尾添加 const
        return period == other.period && initialInventory == other.initialInventory;
    }

    // 允许哈希结构体访问私有成员
    // friend struct
    friend struct std::hash<State>;

    // define operator < or give a self defined comparator for sorting map
    bool operator<(const State &other) const {
        if (period < other.period) {
            return true;
        }
        if (period == other.period) {
            if (initialInventory < other.initialInventory) {
                return true;
            }
            return false;
        }
        return false;
    }
};


// `std::hash<State>` 需要特化
template<> // 表示模版特化， override 标准库中的 hash 生成函数
struct std::hash<State> {
    // size_t 表示无符号整数
    size_t operator()(const State &s) const noexcept {
        // noexcept 表示这个函数不会抛出异常
        // boost 的哈希计算更安全
        std::size_t seed = 0;
        boost::hash_combine(seed, s.period);
        boost::hash_combine(seed, s.initialInventory);
        return seed;

        // return std::hash<int>()(s.period) ^ std::hash<double>()(s.initialInventory) << 1; // 计算哈希值
        // std::hash<int>() 是一个 std::hash<int> 类型的对象，调用 () 运算符可以计算 obj.id（整数）的哈希值
        // ^（异或）是位运算，不会造成进位，适合合并多个哈希值
        // 这里的 << 1 左移 1 位（相当于乘 2），让哈希值更加分散，避免简单叠加导致哈希冲突
    }
};

State::State() = default;

State::State(const int period, const double initialInventory): period(period), initialInventory(initialInventory) {
};

double State::getInitialInventory() const {
    return initialInventory;
}

int State::getPeriod() const {
    return period;
}

void State::print() const {
    std::cout << "period: " << period << ", ini I: " << initialInventory << std::endl;
}


class ProbabilityMassFunctions {
    double truncatedQuantile;
    double stepSize;
    std::string distributionName;

public:
    ProbabilityMassFunctions(double truncatedQuantile, double stepSize, std::string distributionName);

    // std::string getName();

    void checkName() const;

    static double poissonPMF(int k, double lambda);

    [[nodiscard]] std::vector<std::vector<std::vector<double> > > getPMF(std::span<double> demands) const;

    [[nodiscard]] std::vector<std::vector<std::vector<double> > >
    getPMFPoisson(std::span<double> demands) const;

    static int poissonQuantile(double p, double lambda);

    static double poissonCDF(int k, double lambda);
};

// initializing the class
ProbabilityMassFunctions::ProbabilityMassFunctions(
    const double truncatedQuantile, const double stepSize, std::string distributionName)
    : truncatedQuantile(truncatedQuantile), stepSize(stepSize), distributionName(std::move(distributionName)) {
    checkName();
} // std::move for efficiency passing in string and vector

void ProbabilityMassFunctions::checkName() const {
    auto name = distributionName;
    std::ranges::transform(name, name.begin(), ::tolower);
    if (name != "poisson") {
        std::cout << " distribution not found or to do next for this distribution\n";
        raise(-1);
    }
}

// get the probability mass function value of Poisson
double ProbabilityMassFunctions::poissonPMF(const int k, const double lambda) {
    if (k < 0 || lambda <= 0) return 0.0; // 确保参数合法
    return (std::pow(lambda, k) * std::exp(-lambda)) / std::tgamma(k + 1);
    // tgamma(k+1) is a gamma function, 等同于factorial(k)
}


// get cumulative distribution function value of Poisson
double ProbabilityMassFunctions::poissonCDF(const int k, const double lambda) {
    double cumulative = 0.0;
    double term = std::exp(-lambda);
    for (int i = 0; i <= k; ++i) {
        cumulative += term;
        if (i < k)
            term *= lambda / (i + 1); // 递推计算 P(X=i)
    }

    return cumulative;
}

// get inverse cumulative distribution function value of Poisson
int ProbabilityMassFunctions::poissonQuantile(const double p, const double lambda) {
    int low = 0, high = std::max(100, static_cast<int>(lambda * 3)); // 初始搜索区间
    while (low < high) {
        if (const int mid = (low + high) / 2; poissonCDF(mid, lambda) < p) {
            low = mid + 1;
        } else {
            high = mid;
        }
    }
    return low;
}

// get probability mass function values for each period of Poisson
std::vector<std::vector<std::vector<double> > > ProbabilityMassFunctions::
getPMF(const std::span<double> demands) const {
    if (distributionName == "poisson") {
        return getPMFPoisson(demands);
    }
    return {};
}

// get probability mass function values for each period of Poisson
std::vector<std::vector<std::vector<double> > > ProbabilityMassFunctions::
getPMFPoisson(const std::span<double> demands) const {
    const auto T = demands.size();
    int supportLB[T];
    int supportUB[T];
    for (int i = 0; i < T; ++i) {
        supportUB[i] = poissonQuantile(truncatedQuantile, demands[i]);
        supportLB[i] = poissonQuantile(1 - truncatedQuantile, demands[i]);
    }
    std::vector<std::vector<std::vector<double> > > pmf(T, std::vector<std::vector<double> >());
    for (int t = 0; t < T; ++t) {
        const int demandLength = static_cast<int>((supportUB[t] - supportLB[t] + 1) / stepSize);
        pmf[t] = std::vector<std::vector<double> >(demandLength, std::vector<double>());
        for (int j = 0; j < demandLength; ++j) {
            pmf[t][j] = std::vector<double>(2);
            pmf[t][j][0] = supportLB[t] + j * stepSize;
            const int demand = static_cast<int>(pmf[t][j][0]);
            pmf[t][j][1] = poissonPMF(demand, demands[t]) / (2 * truncatedQuantile - 1);
        }
    }

    return pmf;
}

class NewsvendorDP {
    int T;
    int capacity;
    double stepSize;
    double fixOrderCost;
    double unitVariOrderCost;
    double unitHoldCost;
    double unitPenaltyCost;
    double truncatedQuantile;
    double max_I;
    double min_I;

    std::vector<std::vector<std::vector<double> > > pmf;

     std::unordered_map<State, double> cacheActions{};
     std::unordered_map<State, double> cacheValues{};

//    std::map<State, double> cacheActions{};
//    std::map<State, double> cacheValues{};

public:
    NewsvendorDP(size_t T, int capacity, double stepSize, double fixOrderCost, double unitVariOrderCost,
                 double unitHoldCost, double unitPenaltyCost, double truncatedQuantile, double max_I, double min_I,
                 std::vector<std::vector<std::vector<double> > > pmf);

    [[nodiscard]] std::vector<double> feasibleActions() const;

    [[nodiscard]] State stateTransitionFunction(const State &state, double action, double demand) const;

    [[nodiscard]] double immediateValueFunction(const State &state, double action, double demand) const;

    [[nodiscard]] double getOptAction(const State &tate);

    [[nodiscard]] auto getTable() const;

    double recursion(const State &state);
};

NewsvendorDP::NewsvendorDP(const size_t T, const int capacity,
                           const double stepSize, const double fixOrderCost,
                           const double unitVariOrderCost,
                           const double unitHoldCost, const double unitPenaltyCost,
                           const double truncatedQuantile, const double max_I,
                           const double min_I,
                           std::vector<std::vector<std::vector<double> > > pmf): T(static_cast<int>(T)),
    capacity(capacity),
    stepSize(stepSize),
    fixOrderCost(fixOrderCost),
    unitVariOrderCost(unitVariOrderCost),
    unitHoldCost(unitHoldCost), unitPenaltyCost(unitPenaltyCost), truncatedQuantile(truncatedQuantile),
    max_I(max_I), min_I(min_I), pmf(std::move(pmf)) {
};


std::vector<double> NewsvendorDP::feasibleActions() const {
    const int QNum = static_cast<int>(capacity / stepSize);
    std::vector<double> actions(QNum);
    for (int i = 0; i < QNum; i = i + 1) {
        actions[i] = i * stepSize;
    }
    return actions;
}

State NewsvendorDP::stateTransitionFunction(const State &state, const double action, const double demand) const {
    double nextInventory = state.getInitialInventory() + action - demand;
    if (state.getPeriod() == 1) {
        (void) nextInventory;
    }
    if (nextInventory > 0) {
        (void) nextInventory;
    }
    nextInventory = nextInventory > max_I ? max_I : nextInventory;
    nextInventory = nextInventory < min_I ? min_I : nextInventory;

    const int nextPeriod = state.getPeriod() + 1;
    // C++11 引入了统一的列表初始化（Uniform Initialization），鼓励使用大括号 {} 初始化类
    const auto newState = State{nextPeriod, nextInventory};

    return newState;
}

double NewsvendorDP::immediateValueFunction(const State &state, const double action, const double demand) const {
    const double fixCost = action > 0 ? fixOrderCost : 0;
    const double variCost = action * unitVariOrderCost;
    double nextInventory = state.getInitialInventory() + action - demand;
    nextInventory = nextInventory > max_I ? max_I : nextInventory;
    nextInventory = nextInventory < min_I ? min_I : nextInventory;
    const double holdCost = std::max(unitHoldCost * nextInventory, 0.0);
    const double penaltyCost = std::max(-unitPenaltyCost * nextInventory, 0.0);

    const double totalCost = fixCost + variCost + holdCost + penaltyCost;
    return totalCost;
}

double NewsvendorDP::getOptAction(const State &state) {
    return cacheActions[state];
}

auto NewsvendorDP::getTable() const {
    size_t stateNums = cacheActions.size();
    std::vector<std::vector<double> > table(stateNums, std::vector<double>(3));
    int index = 0;
    for (const auto &[fst, snd]: cacheActions) {
        table[index][0] = fst.getPeriod();
        table[index][1] = fst.getInitialInventory();
        table[index][2] = snd;
        index++;
    }
    return table;
}

double NewsvendorDP::recursion(const State &state) {
    double bestQ = 0.0;
    double bestValue = std::numeric_limits<double>::max();
    const std::vector<double> actions = feasibleActions();
    for (const double action: feasibleActions()) {
        double thisValue = 0;
        for (auto demandAndProb: pmf[state.getPeriod() - 1]) {
            thisValue += demandAndProb[1] * immediateValueFunction(state, action, demandAndProb[0]);
            if (state.getPeriod() < T) {
                auto newState = stateTransitionFunction(state, action, demandAndProb[0]);
                (void) action;
                if (cacheValues.contains(newState)) {
                    // some issues here
                    thisValue += demandAndProb[1] * cacheValues[newState];
                } else {
                    thisValue += demandAndProb[1] * recursion(newState);
                }
            }
        }
        if (thisValue < bestValue) {
            bestValue = thisValue;
            bestQ = action;
        }
    }
    cacheActions[state] = bestQ;
    cacheValues[state] = bestValue;
    return bestValue;
}


int main() {
    std::vector<double> demands(30, 20);
    const std::string distribution_type = "poisson";
    constexpr int capacity = 100; // maximum ordering quantity
    constexpr double stepSize = 1.0;
    constexpr double fixOrderCost = 0;
    constexpr double unitVariOderCost = 1;
    constexpr double unitHoldCost = 2;
    constexpr double unitPenaltyCost = 10;
    constexpr double truncQuantile = 0.9999; // truncated quantile for the demand distribution
    constexpr double maxI = 500; // maximum possible inventory
    constexpr double minI = -300; // minimum possible inventory


    const auto pmf = ProbabilityMassFunctions(truncQuantile, stepSize, distribution_type).getPMF(demands);
    const size_t T = demands.size();
    auto model = NewsvendorDP(T, capacity, stepSize, fixOrderCost, unitVariOderCost, unitHoldCost, unitPenaltyCost,
                              truncQuantile, maxI, minI, pmf);

    const auto initialState = State(1, 0);
    const auto start_time = std::chrono::high_resolution_clock::now();
    const auto optValue = model.recursion(initialState);
    const auto end_time = std::chrono::high_resolution_clock::now();
    const std::chrono::duration<double> duration = end_time - start_time;
    std::cout << "planning horizon is " << T << " periods" << std::endl;
    std::cout << "running time of C++ is " << duration << std::endl;
    std::cout << "Final optimal value is: " << optValue << std::endl;
     const auto optQ = model.getOptAction(initialState);
     std::cout << "Optimal Q is: " << optQ << std::endl;
    // auto table = model.getTable();
    return 0;
}

Java 代码

import java.util.*;
import java.util.concurrent.ConcurrentHashMap;
import java.util.concurrent.ConcurrentSkipListMap;
import java.util.function.Function;
import java.util.stream.DoubleStream;
import java.util.stream.IntStream;

public class CLSP {
    double[][][] pmf;

    public CLSP(double[][][] pmf) {
        this.pmf = pmf;
    }

    class State {
        int period;
        double initialInventory;

        public State(int period, double initialInventory) {
            this.period = period;
            this.initialInventory = initialInventory;
        }

        public double[] getFeasibleActions() {
            return actionGenerator.apply(this);
        }

        @Override
        public int hashCode() {
            String hash = "";
            hash = hash + period + initialInventory;
            return hash.hashCode();
        }

        @Override
        public boolean equals(Object o) {
            if (o instanceof State)
                return ((State) o).period == this.period &&
                        ((State) o).initialInventory == this.initialInventory;
            else
                return false;
        }

        @Override
        public String toString() {
            return "period = " + period + ", " + "initialInventory = " + initialInventory;
        }
    }

    Function<State, double[]> actionGenerator;


    interface StateTransitionFunction<S, A, R, S2> {
        public S2 apply(S s, A a, R r);
    }

    StateTransitionFunction<State, Double, Double, State> stateTransition;


    interface ImmediateValueFunction<S, A, R, V> {
        public V apply(S s, A a, R r);
    }

    ImmediateValueFunction<State, Double, Double, Double> immediateValue;

    Comparator<State> keyComparator = (o1, o2) -> o1.period > o2.period ? 1 :
            o1.period == o2.period ? Double.compare(o1.initialInventory, o2.initialInventory) : -1;

//
    ConcurrentSkipListMap<State, Double> cacheActions = new ConcurrentSkipListMap<>(keyComparator);
    ConcurrentSkipListMap<State, Double> cacheValues = new ConcurrentSkipListMap<>(keyComparator);

    double f(State state) {
        return cacheValues.computeIfAbsent(state, s -> {
//            double val = Arrays.stream(s.getFeasibleActions())
//                    .map(orderQty -> Arrays.stream(pmf[s.period - 1])
//                            .mapToDouble(p -> p[1] * immediateValue.apply(s, orderQty, p[0]) +
//                                    (s.period < pmf.length ?
//                                            p[1] * f(stateTransition.apply(s, orderQty, p[0])) : 0))
//                            .sum())
//                    .min()
//                    .getAsDouble();
//            double bestOrderQty = Arrays.stream(s.getFeasibleActions())
//                    .filter(orderQty -> Arrays.stream(pmf[s.period - 1])
//                            .mapToDouble(p -> p[1] * immediateValue.apply(s, orderQty, p[0]) +
//                                    (s.period < pmf.length ?
//                                            p[1] * f(stateTransition.apply(s, orderQty, p[0])) : 0))
//                            .sum() == val)
//                    .findAny()
//                    .getAsDouble();
//            cacheActions.putIfAbsent(s, bestOrderQty);
//            return val;
//        });
//    }

            double[] feasibleActions = state.getFeasibleActions();
            double[][] dAndP = pmf[state.period - 1]; // demandAndPossibility
            double[] QValues = new double[feasibleActions.length];
            double val = Double.MAX_VALUE;

            double bestOrderQty = 0;
            for (int i = 0; i < feasibleActions.length; i++) {
                double orderQty = feasibleActions[i];
                double thisQValue = 0;
                for (int j = 0; j < dAndP.length; j++) {
                    thisQValue += dAndP[j][1] * immediateValue.apply(state, orderQty, dAndP[j][0]);
                    if (state.period < pmf.length) {
                        State newState = stateTransition.apply(state, orderQty, dAndP[j][0]);
                        thisQValue += dAndP[j][1] * f(newState);
                        }
                    }
                QValues[i] = thisQValue;

                if (QValues[i] < val) {
                    val = QValues[i];
                    bestOrderQty = orderQty;
                }

            }
            this.cacheActions.putIfAbsent(state, bestOrderQty);
//            cacheValues.put(state, val);
            return val;
        });
    }



    public static void main(String[] args) {


        double initialInventory = 0;
        double[] meanDemand = new double[30];
        Arrays.fill(meanDemand, 20);

        double truncationQuantile = 0.9999;
        double stepSize = 1;
        double minState = -150;
        double maxState = 300;
        int T = meanDemand.length;

        double fixedOrderingCost = 0;
        double proportionalOrderingCost = 1;
        double holdingCost = 2;
        double penaltyCost = 10;

        int maxOrderQuantity = 100;

        Distribution[] distributions = IntStream.iterate(0, i -> i + 1)
                .limit(T)
                .mapToObj(i -> new PoissonDist(meanDemand[i]))
//	                                              .mapToObj(i -> new UniformDist(0, meanDemand[i]))
                //.mapToObj(i -> new NormalDist(meanDemand[i], 0.25 * meanDemand[i]))
                .toArray(Distribution[]::new); // replace for loop
        double[] supportLB = IntStream.iterate(0, i -> i + 1)
                .limit(T)
                .mapToDouble(i -> distributions[i].inverseF(1 - truncationQuantile))
                .toArray();
        double[] supportUB = IntStream.iterate(0, i -> i + 1)
                .limit(T)
                .mapToDouble(i -> distributions[i].inverseF(truncationQuantile))
                .toArray();
        double[][][] pmf = new double[T][][];
        for (int i = 0; i < T; i++) {
            int demandLength = (int) ((supportUB[i] - supportLB[i] + 1) / stepSize);
            pmf[i] = new double[demandLength][];
            // demand values are all integers
            for (int j = 0; j < demandLength; j++) {
                pmf[i][j] = new double[2];
                pmf[i][j][0] = supportLB[i] + j * stepSize;
                int demand = (int) pmf[i][j][0];
                if (distributions[0] instanceof DiscreteDistribution) {
                     // double probabilitySum = distributions[i].cdf(supportUB[i]) - distributions[i].cdf(supportLB[i]);
                    double probabilitySum = 2 * truncationQuantile - 1;
                    pmf[i][j][1] = ((DiscreteDistribution) distributions[i]).prob(demand) / probabilitySum;
                } else {
                    double probabilitySum = distributions[i].cdf(supportUB[i] + 0.5 * stepSize)
                            - distributions[i].cdf(supportLB[i] - 0.5 * stepSize);
                    pmf[i][j][1] = (distributions[i].cdf(pmf[i][j][0] + 0.5 * stepSize)
                            - distributions[i].cdf(pmf[i][j][0] - 0.5 * stepSize)) / probabilitySum;
                }
            }
        }

        CLSP inventory = new CLSP(pmf);

        inventory.actionGenerator = s -> {
            return DoubleStream.iterate(0, i -> i + stepSize).limit(maxOrderQuantity + 1).toArray();
        };

        inventory.stateTransition = (state, action, randomDemand) -> {
            double nextInventory = state.initialInventory + action - randomDemand;
            nextInventory = nextInventory > maxState ? maxState : nextInventory;
            nextInventory = nextInventory < minState ? minState : nextInventory;
            return inventory.new State(state.period + 1, nextInventory);
        };


        inventory.immediateValue = (state, action, randomDemand) ->
        {
            double fixedCost = action > 0 ? fixedOrderingCost : 0;
            double variableCost = proportionalOrderingCost * action;
            double inventoryLevel = state.initialInventory + action - randomDemand;
            double holdingCosts = holdingCost * Math.max(inventoryLevel, 0);
            double penaltyCosts = penaltyCost * Math.max(-inventoryLevel, 0);
            double totalCosts = fixedCost + variableCost + holdingCosts + penaltyCosts;
            return totalCosts;
        };

        int period = 1;
        State initialState = inventory.new State(period, initialInventory);
        long currTime2 = System.currentTimeMillis();

        double finalValue = inventory.f(initialState);
        double time = (System.currentTimeMillis() - currTime2) / 1000.000;
        System.out.println("planning horizon is " + meanDemand.length + " periods");
        System.out.println("running time of Java is " + time + " s");

        System.out.println("final optimal expected value is: " + finalValue);

        double optQ = inventory.cacheActions.get(inventory.new State(period, initialInventory));
        System.out.println("optimal order quantity in the first priod is : " + optQ);

    }
}