一种核电阀门寿命试验中的最佳应力和样本分配确定方法
技术领域
本发明属于核电阀门可靠性分析领域，尤其争对阀门加速寿命试验。
背景技术
传统的核电阀门加速寿命试验存在着试验时间较长、样本量较大的缺点，因此，有必要对加速寿命试验方案进行优化设计，寻找一组最佳的应力水平和试样分配比例，使阀门在正常应力水平时寿命估计值的精度最高，试验时间缩短，同时减少样本数量。目前常见的全局最优求解法有拉格朗日法、线性规划法、以及一些人工智能算法比如遗传算法、粒子群算法、模拟退火算法等。优化算法已经在科技、机械工程领域十分普及，但是往往现实中的优化问题都伴随着各种约束条件，因此衍生出了基于遗传算法的带约束求解最优法等，由于传统惩罚函数中惩罚系数的范围很难确定，因此优化难度就加大了。
发明内容
本发明对现有核电阀门加速寿命试验方案设计领域技术背景的不足，根据试验的加速性和试验数据的统计精度要求，确定试验优化的约束条件，本发明采用内点罚函数法对试验优化的约束问题进行处理。
本发明技术方案为一种核电阀门加速寿命试验中的最佳应力和样本分配确定方法，该方法包含如下：
步骤1：选定应用于核电阀门的具体应力类型为温度应力，确定阀门寿命分布服从威布尔分布，且威布尔分布形状参数为m，尺度参数为η；于是寿命分布参数σ＝1/m，μ＝lnη，确定该应力x与核电阀门寿命分布参数μ的线性关系μ＝γ0+γ1x，并且特征寿命与试验温度之间满足阿伦尼斯方程η＝exp(γ0+γ1x)，γ0，γ1和σ为根据先验数据得到的模型参数；
步骤2：建立优化目标函数，具体步骤如下：
步骤2.1：在核电阀门加速寿命试验方案优化设计中通常以线性-极值模型下应力水平个数为M时的方差因子VM为优化目标函数，求解最佳应力水平和样本分配比。
步骤2.2：建立目标函数为：
s.t.g(i)＝Ti-1-Ti＜0,(i＝1,2,...,M)
其中，VM表示应力水平为M时的方差因子，x0表示正常应力水平，n表示总样本个数，zp表示标准极值分布的P阶分位数，g(i)表示温度应力约束，h(i)表示样本约束，Ti表示温度应力，pi表示样本分配比，M表示应力水平个数，F表示方差因子中关于应力水平xi和样本分配比pi的信息矩阵；
步骤3：针对核电阀门加速寿命试验方案优化设计问题，采用内点罚函数法对试验优化的约束问题进行处理，求解目标函数得到最佳温度应力水平xi和样本分配比pi。
进一步的，步骤3的具体方法为：
步骤3.1：对优化目标函数引入惩罚函数，将不等式约束函数：
s.t.g(i)＝Ti-1-Ti＜0,(i＝1,2,...,M)
转化为新目标函数：
其中，r(k)表示取惩罚因子，r(k+1)＝C·r(k)，C为惩罚因子缩小系数；
步骤3.2：计算目标函数最小值时的应力水平和最佳样本分配个数。
进一步的，所述步骤3.2的计算方法为：
步骤3.2.1：确定目标函数f(x)，假设目标函数具有多个决策变量，即x＝(x1,x2,...,xn)；
步骤3.2.2：优化前确定迭代次数N，初始步长λ，优化变量的初始值x0以及控制精度ξ，其中ξ＞0；
步骤3.2.3：判断是否达到得设定的迭代次数，即k＜N是否满足，k＝1。
步骤3.2.4：若k＝1进行第一次迭代，则在(-σ,σ)之间随机生成n个n维向量μi＝(μi1,μi2,...,μin)，其中σ＝xmax-xmin，然后将μi进行转化得到令xi＝x+σμ′i，得到{x1,x2,...,xn}，代入目标函数值进行计算，得到使目标函数值最小的一组解，完成第一次的迭代；
步骤3.2.5：将带约束的目标函数转化为无约束的目标函数，然后对目标函数进行计算，如果满足条件f1(x1)＜f1(x)，那么就找到了比起始值好的点，使得k＝1，将x＝x1，返回步骤3.2.3，若不满足条件，则k＝k+1，r(k+1)＝C·r(k)返回步骤3.2.4；
步骤3.2.6：如果迭代次数已经满足条件，仍然未找到最优的值，那么输出当前的解，但并非最优解；此时达到终止条件的时候，即σ＜ξ，算法截至；不满足截止条件则继续迭代。
本发明有益效果为：通过试验分析，在进行核电阀门加速寿命试验优化设计时将随机游走优化算法和惩罚函数法结合有效的解决了带约束的优化求解问题，寻找到一组最佳的应力水平和试样分配比例，并且在求变量更优值的同时使得惩罚因子的确定更加的准确和简易，为其他优化算法求解带约束问题提供了新思路。
说明书附图
图1为随机游走带约束优化求解流程图
图2为改进随机游走优化算法流程图
具体实施方式
步骤1：选定应用于核电阀门的具体应力类型为温度应力，确定阀门寿命分布服从威布尔分布，且威布尔分布形状参数为m，尺度参数为η。于是寿命分布参数σ＝1/m，μ＝lnη，确定该应力x与核电阀门寿命分布参数μ的线性关系μ＝γ0+γ1x，并且特征寿命与试验温度之间满足阿伦尼斯方程η＝exp(γ0+γ1x)，γ0，γ1和σ为根据先验数据得到的模型参数；
步骤2：建立优化目标函数，具体步骤如下：
步骤2.1：在核电阀门加速寿命试验方案优化设计中通常以应力水平为M时的方差因子VM为优化目标函数，求解最佳应力水平和样本分配比。
步骤2.2：建立目标函数为：
s.t.g(i)＝Ti-1-Ti＜0,(i＝1,2,...,M)
其中，VM表示应力水平为M时的方差因子，x0表示正常应力水平，n表示总样本个数，zp表示标准极值分布的P阶分位数，g(i)表示温度应力约束，h(i)表示样本约束，Ti表示温度应力，pi表示样本分配比，M表示应力水平个数，F表示方差因子中关于应力水平xi和样本分配比pi的信息矩阵；
步骤3：针对核电阀门加速寿命试验方案优化设计问题，采用内点罚函数法对试验优化的约束问题进行处理，求解目标函数得到最佳应力水平xi和样本分配比pi。
进一步的，步骤3的具体方法为：
步骤3.1：对优化目标函数引入惩罚函数，将不等式约束函数：
s.t.g(i)＝Ti-1-Ti＜0,(i＝1,2,...,M)
转化为新目标函数：
其中，r(k)表示取惩罚因子，r(k+1)＝C·r(k)，C为惩罚因子缩小系数；
步骤3.2：计算目标函数最小值时的应力水平和最佳样本分配个数。
进一步的，所述步骤3.2的计算方法为：
步骤3.2.1：确定目标函数f(x)，假设目标函数具有多个决策变量，即x＝(x1,x2,...,xn)；
步骤3.2.2：优化前确定迭代次数N，初始步长λ，优化变量的初始值x0以及控制精度ξ，其中ξ＞0；
步骤3.2.3：判断是否达到得设定的迭代次数，即k＜N是否满足，k＝1。
步骤3.2.4：若k＝1进行第一次迭代，则在(-σ,σ)之间随机生成n个n维向量μi＝(μi1,μi2,...,μin)，其中σ＝xmax-xmin，然后将μi进行转化得到令xi＝x+σμ′i，得到{x1,x2,...,xn}，代入目标函数值进行计算，得到使目标函数值最小的一组解，完成第一次的迭代；
步骤3.2.5：将带约束的目标函数转化为无约束的目标函数，然后对目标函数进行计算，如果满足条件f1(x1)＜f1(x)，那么就找到了比起始值好的点，使得k＝1，将x＝x1，返回步骤3.2.3，若不满足条件，则k＝k+1，r(k+1)＝C·r(k)返回步骤3.2.4；
步骤3.2.6：如果迭代次数已经满足条件，仍然未找到最优的值，那么输出当前的解，但并非最优解；此时达到终止条件的时候，即σ＜ξ，算法截至；不满足截止条件则继续迭代。
本发明提出的优化算法计算步骤如下：
步骤1：在阀门加速寿命试验优化设计中，假设加速试验有2温度水平T1＜T2，确定试验中样品的分配比例和极限应力水平，满足p1+p2＝1。VM表示应力水平为M时的方差因子，g(i)和h(i)为试验变量实际取值范围约束。
目标函数为：
s.t.g(i)＝Ti-1-Ti＜0,(i＝1,2)
步骤2：根据摸底试验可以得到模型参数的估计值zp＝0.5；根据试验条件可以知道，该产品的正常工作应力水平T0＝85℃，其转化应力为x0＝1000/(273.15+T0)，极限工作应力水平，Tmax＝140℃，其转化应力为给定置信度为γ＝60％，则Kγ＝0.84162，然后置信区间宽度为2W＝0.4，根据式子n＝VM(Kγσ/W)2求出样本量，引入惩罚函数得到新目标函数为：
式子中以温度应力水平和样本比例为决策变量，x：T1,T2,p1,p2；r(1)＞r(2)＞r(3)＞...＞r(k)为惩罚因子，一般取正数，
为了构造递减的惩罚因子序列，则
r(k+1)＝Cr(k)
式子中，C∈(0,1)为惩罚因子缩小系数。一般情况下，取惩罚因子r(1)＝1，惩罚因子缩小系数C＝0.5～0.7。
步骤3：计算目标函数最小值时的温度应力水平值和最佳样本分配比例。
进一步的，所述步骤3的计算方法为：
步骤3.1.1：确定目标函数f(x)，假设目标函数具有多个决策变量，即x：T1,T2,p1,p2；
步骤3.1.2：优化前确定迭代次数N＝1000，初始步长λ＝0.5，优化变量的初始值x0:[2.793,2.421,0,1,1]以及控制精度ξ＝0.00001(其中ξ＞0，一般取值很小)；
步骤3.1.3：判断是否达到得设定的迭代次数，即k＜N是否满足，k＝1。
步骤3.1.4：若k＝1进行第一次迭代，则在([-0.372,-1],[0.372,1])之间随机生成10个向量μi＝(μi1,μi2,...,μin)，然后将μi进行转化得到令xi＝x+σμ′i，得到{x1,x2,...,x10}，代入目标函数值进行计算，得到使目标函数值最小的一组解，完成第一次的迭代。
步骤3.1.5：将带约束的目标函数转化为无约束的目标函数，然后对目标函数进行计算，如果满足条件f1(x1)＜f1(x)，那么就找到了比起始值好的点，使得k＝1，将x＝x1，返回步骤3.1.3，若不满足条件，则k＝k+1，r(k+1)＝Cr(k)，C＝0.5～0.7，返回步骤3.1.4；
步骤3.1.6：如果迭代次数已经满足条件，仍然未找到最优的值，那么输出当前的解，但并非最优解。此时达到终止条件的时候，即σ＜ξ，算法截至；不满足截止条件则继续迭代。
表1传统的随机游走和改进随机游走结果对比
