开发者选择生成对抗网络（GAN）学习固体火箭的内弹道数据，目的是尝试让AI生成固体火箭的内弹道数据以供设计参考使用。

数据由warmonkey于2011年3月开源的仿真程序“burnsim”生成：

“burnsim”程序的原作者声明：（程序）已经经过试车检验，可认为计算较准确。

开发者认为：AI设计像火箭一类的工业产品面临的一大问题便是高品质训练数据短缺，在当下（2023年），文本数据和图像数据的获取“几乎没有成本”，因此相关领域在应用机器学习技术之后得到飞速发展。至于工业产品设计，获取数据就很困难了——以被爱好者“做烂”的固体火箭为例，试问哪一个爱好者能够通过实际点火试车获得大量可靠的内弹道数据？在数据获取成本高的情况下，只能考虑使用可靠的仿真程序生成训练使用的数据。

Warmonkey给出的仿真程序是m程序，包含3个m文件：

机器学习一般使用py程序，为此，开发者使用“XXXXXXXXXXXXXXXp”网站部署的人工智能，将m程序转化为py程序，如图：

这个网站转换程序之后不需要太多修正就能直接运行（py程序依赖的函数库需要手动补全），将绘制图像的代码去除，并对程序进行适当的修正，运行XXXXXXXXXX：

# BurnSim release 1.3
# 固体火箭发动机内弹道仿真器
import math
import numpy as np
from openpyxl import load_workbook

def burnsim_propellant(pressure, propname, argcname):
    if propname == 'KNSU':
        if argcname == 'vc':
            result = 923  # m/s 特征速度
        elif argcname == 'density':
            result = 1.81  # g/cm^3 密度
        elif argcname == 'burnspeed':
            a, n = 8.263, 0.319  # 燃速系数
            result = a * pressure ** n  # 燃速 r = a * pressure ^ n 燃速指数
        elif argcname == 'k':
            result = 1.044  # 比热比
    elif propname == 'KNDX':
        if argcname == 'vc':  # m/s 特征速度
            result = 897  # 891,916
        elif argcname == 'density':
            result = 1.879 * 0.94  # g/cm^3 密度
        elif argcname == 'burnspeed':
            if pressure < 0.779:  # 0.1MPa以下按此估算
                a, n = 8.875, 0.619  # 燃速系数
            elif pressure < 2.572:
                a, n = 7.553, -0.009  # 燃速系数
            elif pressure < 5.930:
                a, n = 3.841, 0.688  # 燃速系数
            elif pressure < 8.502:
                a, n = 17.20, -0.148  # 燃速系数
            else:  # 11.50MPa以上按此估算
                a, n = 4.775, 0.442  # 燃速系数
            result = a * pressure ** n  # 燃速
        elif argcname == 'k':
            result = 1.043  # 比热比
    elif propname == 'KNSB':
        if argcname == 'vc':  # m/s 特征速度
            result = 914
        elif argcname == 'density':
            result = 1.78  # g/cm^3 密度
        elif argcname == 'burnspeed':
            if pressure < 0.807:  # 0.1MPa以下按此估算
                a, n = 10.708, 0.625  # 燃速系数
            elif pressure < 1.503:
                a, n = 8.763, -0.314  # 燃速系数
            elif pressure < 3.792:
                a, n = 7.852, -0.013  # 燃速系数
            elif pressure < 7.033:
                a, n = 3.907, 0.535  # 燃速系数
            else:  # 10.67MPa以上按此估算
                a, n = 9.653, 0.064  # 燃速系数
            result = a * pressure ** n  # 燃速
        elif argcname == 'k':
            result = 1.042  # 比热比
    # 此处添加新燃料
    return result
def burnsim(prop,d_nozzle_begin,D_outer,d_hole,L):
    solution_name='dataset'   #设计方案名称
#    d_nozzle_begin=20  #喷喉初始直径，单位毫米
    nozzle_burn_speed=0   #喷管烧蚀速率
#    args = parser.parse_args()


#    d_hole = [8]  # mm 药柱内孔径
#    D_outer = 45  # mm 药柱外径
    N = 5  # 药柱能够燃烧的横断面数
#    L = [300]  # mm 药柱可燃部分长度
#    prop = 'KNDX'  # 推进剂名称
    ap = 0.1  # MPa 大气压
    n = 0.85  # 喷管效率
    g = 9.8  # m/s^2 重力加速度

    # 初始值 单位
    pi = math.pi  # 圆周率
    d_nozzle = [d_nozzle_begin]  # mm
    r = [burnsim_propellant(ap, prop, 'burnspeed')]  # mm/s
    V_gram = [(D_outer ** 2 - d_hole[0] ** 2) * pi / 4 * L[0] / 1000]  # cm^3
    c = burnsim_propellant(0, prop, 'vc')  # m/s
    de = burnsim_propellant(0, prop, 'density')  # g/cm^3
    k = burnsim_propellant(0, prop, 'k')  # 比热比
    I = 0  # N*s
    t = 0  # s
    dt = 0.001  # s 分析时间片长度,数值小精度高,数值大需要内存少,修改后注意要更改图像的X坐标名称
    
    i = 1  # 循环计数
    while V_gram[i - 1] > 0:
        # 计算公式 得数含义 单位
        Ae = d_nozzle[i - 1] ** 2 * pi / 4  # 当前喷管直径 mm^2
        S_burn = (D_outer ** 2 - d_hole[i - 1] ** 2) * pi / 4 * N + d_hole[i - 1] * pi * L[i - 1]  # 当前燃烧面积 mm^2
        Kn = S_burn / Ae  # 喷燃比
        P = Kn * r[i - 1] * c * de / 1e6  # 压力 MPa
        r.append(burnsim_propellant(P, prop, 'burnspeed'))  # 当前燃速 mm/s
        k_thrust = n * abs((2 * k ** 2 / (k - 1) * (2 / (k + 1)) ** ((k + 1) / (k - 1))) ** 0.5 * (1 - (ap / P) ** ((k - 1) / k)))  # 当前推力系数
        F = k_thrust * Ae * P  # 推力
        I += F * dt  # 总冲
        d_hole.append(d_hole[i - 1] + 2 * r[i] * dt)  # 下一时刻的药柱内孔径 mm
        L.append(L[i - 1] - N * dt * r[i])  # 下一时刻的药柱长度 mm
        V_gram.append((D_outer ** 2 - d_hole[i] ** 2) * pi / 4 * L[i] / 1000)  # 下一时刻的药柱体积 cm^3
        d_nozzle.append(d_nozzle[i - 1] + 2 * nozzle_burn_speed * dt)  # 下一时刻的喷管内径 mm
    
        i += 1
        t += dt

    Isp = I / 9.8 / V_gram[0] / de * 1000  # 比冲 单位为秒
    F_avg = I / t  # 平均推力
    M_propellant = V_gram[0] * de  # 药柱质量 单位为g
    data=[Ae,S_burn,Kn,P,k_thrust,F,I,Isp,F_avg,M_propellant]
    print(data)
    return data
if __name__ == '__main__':
    prop = 'KNSU'
    d0 = np.linspace(2, 30, 8)
    d1 = np.linspace(35, 50, 8)
    d2 = np.linspace(10, 30, 8)
    d3 = np.linspace(100, 800, 8)
    
    b = []
    
    wb = load_workbook('dataset-KNSU.xlsx')
    sheet = wb.active
    
    # 向excel中写入表头
    sheet['a1'] = '喷管出口面积'
    sheet['b1'] = '燃烧面积'
    sheet['c1'] = '喷燃比'
    sheet['d1'] = '燃烧室压强'
    sheet['e1'] = '推力系数'
    sheet['f1'] = '推力'
    sheet['g1'] = '总冲'
    sheet['h1'] = '比冲'
    sheet['i1'] = '平均推力'
    sheet['j1'] = '药柱质量'
    
    for d00 in d0:
        d_nozzle_begin = d00
        
        for d11 in d1:
            D_outer = d11
            
            for d22 in d2:
                d_hole = [d22]
                
                for d33 in d3:
                    L = [d33]
                    result = burnsim(prop, d_nozzle_begin, D_outer, d_hole, L)
                    b.append(result)
                    
    for bb in b:
        sheet.append(bb)
        
    wb.save('dataset-KNSU.xlsx')
    print("成功写入数据！！！")

自我批评：向excel文件写入生成的数据时只知道嵌套for循环，这样操作增加了程序运行占用的内存空间，不合适。

运行XXXXXXXXXX之后生成一个excel文件“dataset-KNSU.xlsx”；

dataset-KNSU.xlsx 505.80KB XLSX 4次下载

将文件内的数组输入给神经网络展开训练，以下为开发者运行的神经网络程序XXXXXX：

import argparse
import numpy as np
import torchvision.transforms as transforms
from torchvision.utils import save_image
from torch.utils.data import DataLoader,TensorDataset
from torchvision import datasets
from torch.autograd import Variable
import pandas as pd
import torch.nn as nn
import torch.nn.functional as F
import torch
from matplotlib import pyplot as plt

parser = argparse.ArgumentParser()
parser.add_argument("--n_epochs", type=int, default=20, help="number of epochs of training")
parser.add_argument("--batch_size", type=int, default=64, help="size of the batches")
parser.add_argument("--lr", type=float, default=0.02, help="adam: learning rate")
parser.add_argument("--b1", type=float, default=0.5, help="adam: decay of first order momentum of gradient")
parser.add_argument("--b2", type=float, default=0.999, help="adam: decay of first order momentum of gradient")
parser.add_argument("--n_cpu", type=int, default=0, help="number of cpu threads to use during batch generation")
parser.add_argument("--latent_dim", type=int, default=10, help="dimensionality of the latent space")
parser.add_argument("--channels", type=int, default=10, help="number of image channels")
parser.add_argument("--sample_interval", type=int, default=40, help="interval betwen image samples")
opt = parser.parse_args(args=[])
print(opt)

#img_shape = (opt.channels,1,1)

cuda =False
class Generator(nn.Module):
    def __init__(self):
        super(Generator, self).__init__()

        self.model = nn.Sequential(
            nn.Linear(opt.latent_dim, 128),
            nn.LeakyReLU(0.2, inplace=True),
            nn.Linear(128, 256),
            nn.BatchNorm1d(256, 0.8),
            nn.LeakyReLU(0.2, inplace=True),
            nn.Linear(256, 512),
            nn.BatchNorm1d(512, 0.8),
            nn.LeakyReLU(0.2, inplace=True),
            nn.Linear(512, 1024),
            nn.BatchNorm1d(1024, 0.8),
            nn.LeakyReLU(0.2, inplace=True),
            nn.Linear(1024, 1),
            nn.Tanh()
        )

    def forward(self, z):
        img = self.model(z)
        img = img.view(img.size(0), 1)
        return img


class Discriminator(nn.Module):
    def __init__(self):
        super(Discriminator, self).__init__()

        self.model = nn.Sequential(
            nn.Linear(1, 512),
            nn.LeakyReLU(0.2, inplace=True),
            nn.Linear(512, 256),
            nn.LeakyReLU(0.2, inplace=True),
            nn.Linear(256, 10),
            nn.Sigmoid(),
        )

    def forward(self, img):
        img_flat = img.view(img.size(0), -1)
        validity = self.model(img_flat)

        return validity


# Loss function
adversarial_loss = torch.nn.BCELoss()

# Initialize generator and discriminator
generator = Generator()
discriminator = Discriminator()

if cuda:
    generator.cuda()
    discriminator.cuda()
    adversarial_loss.cuda()


# Optimizers
optimizer_G = torch.optim.Adam(generator.parameters(), lr=opt.lr, betas=(opt.b1, opt.b2))
optimizer_D = torch.optim.Adam(discriminator.parameters(), lr=opt.lr, betas=(opt.b1, opt.b2))

Tensor = torch.cuda.FloatTensor if cuda else torch.FloatTensor

# ----------
#  Training
# ----------
dat_loader = pd.read_excel('dataset-KNSU.xlsx', 'Sheet1', skiprows=1, header=None)
dat=dat_loader.values
get_data =torch.from_numpy(dat)
imgs = torch.tensor(get_data)
#print(imgs.shape,type(imgs))
def plot_curve(data):
    fig = plt.figure()
    plt.plot(range(len(data)), data, color='blue')
    plt.legend(['value'], loc='upper right')
    plt.xlabel('step')
    plt.ylabel('value')
    plt.show()
for epoch in range(opt.n_epochs):
    for im in imgs:
        valid = Variable(Tensor(im.size(0), 10).fill_(1.0), requires_grad=False)
        fake = Variable(Tensor(im.size(0), 10).fill_(0.0), requires_grad=False)
        real_imgs = Variable(im.type(Tensor))
        real_imgs=torch.unsqueeze(real_imgs,dim=1)
   #     print(real_imgs.shape)
        optimizer_G.zero_grad()
        z = Variable(Tensor(np.random.normal(0, 1, (im.shape[0], opt.latent_dim))))
#        print(z.shape)
        gen_imgs = generator(z)
#        print(discriminator(gen_imgs).shape)
        g_loss = adversarial_loss(discriminator(gen_imgs), valid)
        g_loss.backward()
       
        optimizer_G.step()
        optimizer_D.zero_grad()
#        nn.utils.clip_grad_norm(generator.parameters, max_norm=1, norm_type=2)
#        nn.utils.clip_grad_norm(discriminator.parameters, max_norm=1, norm_type=2)
        print(discriminator(real_imgs))
        
        real_imgs_normalized = (real_imgs - real_imgs.min()) / (real_imgs.max() - real_imgs.min())
#        print(real_imgs_normalized)
        real_loss = adversarial_loss(discriminator(real_imgs_normalized), valid)

        fake_loss = adversarial_loss(discriminator(gen_imgs.detach()), fake)
        d_loss = (real_loss + fake_loss) / 2
        d_loss.backward()
        optimizer_D.step()
        print(
                "[Epoch %d/%d] [D loss: %f] [G loss: %f]"
                 % (epoch, opt.n_epochs, d_loss.item(), g_loss.item())
              )

运行上面的程序（运行的硬件环境是英特尔i7-7700CPU，8GB内存，不使用CUDA），得到梯度爆炸的结果：

图中在epoch0还没有训练完成的时候，输出的discriminator(real_imgs)出现“NAN”，模型无法收敛，项目暂时停滞。