止于至善

0%

基于粒子群算法的本构参数标定方法

发表于 更新于 分类于

前言

商业有限元软件允许用户定义本构关系，例如ABAQUS可以通过UserMATerial（UMAT）子程序引入复杂的本构关系。复杂本构中的参数往往无法与实验结果建立一一对应的关系，解决参数标定的基本思想是用有限元计算结果去不断逼近实验结果，这会引出两个子问题：

  1. 如何定量描述有限元计算or实验结果之间的误差
  2. 采用什么方法逼近

发扬“愚公移山”精神的试错法（即肉眼衡量+手动迭代）需要用户手动完成有限元的前处理（修改参数）、提交计算、后处理（提取数据）的全部过程，显然十分费时费力。

为了高效地解决这个问题，经调研发现HuLi在其文章中使用了基于粒子群算法来标定参数。个人认为这种方法实现简单、并行计算强、标定效果好，下面介绍我的具体实现方法。

两个子问题

子问题1

前言中的第一个子问题，如何定量描述两种结果之间的误差

在实验中常常得到的是各种载荷-位移曲线，记为\(f_1=g(\delta_1)\)；有限元计算是对实验的模拟，相应地也可以得一条载荷-位移曲线，记为\(f_2=h(\delta_2)\)。事实上，实验和有限元得到数据并不连续，而是众多的离散点，因此连续函数\(~g~\)和\(~h~\)在离散点外的值需要通过插值得到，可以使用MATLAB里的线性插值命令interp1。对于一组材料参数\(\boldsymbol{x}\)，选取\(N\)个数据点，相应的误差可定义为

\[ \tag{1} err(\boldsymbol{x}) = \sqrt{\frac{\sum\nolimits^N_{i=1} \big[g(\delta_i)-h(\delta_i)\big]^2}{N}} \]

子问题2

前言中的第二个子问题，如何逼近，使得误差\(err\)尽量小

如何求式\((1)\)的极小值，在数学上等价为一个最优化问题。粒子群优化算法不需要目标函数的导数信息，且可实现并行计算，故对参数标定这类优化问题有很好的效果。这里简单介绍该算法的基本思想，详细内容可参考blog或《粒子群优化算法》。

粒子群算法认为对每个粒子有两个属性：速度\(\boldsymbol{v}\)和位置\(\boldsymbol{x}\)，其中位置\(\boldsymbol{x}\)在本文的背景下就是一组材料参数。对于每个位置\(\boldsymbol{x}\)，可由式\((1)\)得到\(err\)，并称为该位置的适应度。

在算法开始时，根据指定的速度和位置边界来随机地初始化一系列粒子，即初始粒子群。在第\(~j~\)步迭代过程中，第\(~i~\)个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值\(\mathrm{pFitness}^{(j)}_i\)，并将个体极值与整个粒子群里的其他粒子共享，找到所有迭代步中最优的个体极值作为整个粒子群的当前全局最优解\(\mathrm{gFitness}^{(j)}\)，相应地把个体极值\(\mathrm{pFitness}^{(j)}_i\)对应的位置记为历史最优位置\(\mathrm{pBest}^{(j)}_i\)，把全局最优解\(\mathrm{gFitness}^{(j)}\)对应的位置记为全局最优位置\(\mathrm{gBest}^{(j)}\)。粒子通过跟踪两个最优位置来更新自身，更新公式如下 \[ \tag{2a} \boldsymbol{v}^{(j+1)}_i = \underbrace{\omega^{(j)} \boldsymbol{v}^{(j)}_{i}}_{\text{记忆项}} + \underbrace{c_1~r_1~(\mathrm{pBest}^{(j)}_i - \boldsymbol{x}^{(j)}_i)}_{\text{自身认知项}} + \underbrace{c_2~r_2~(\mathrm{gBest}^{(j)} - \boldsymbol{x}^{(j)}_i)}_{\text{群体认知项}} \]

\[ \tag{2b} \boldsymbol{x}^{(j+1)}_i = \boldsymbol{x}^{(j)}_i + \boldsymbol{v}^{(j+1)}_i \]

式中，\(c_1\)、\(c_2\)为学习因子，一般均取2；\(r_1\)、\(r_2\)为0到1之间的随机数；\(\omega^{(j)}\)为第\(~j~\)次迭代的惯性因子，一般采用常数或线性递减形式。式\((2)\)便是粒子群算法的核心公式，收敛准则可采用\(err(\boldsymbol{x})<e_0\)的简单形式。

程序实现

本文使用ABAQUS有限元软件，采用MATLAB作为主计算平台，负责粒子群算法的实现与有限元算计三个步骤的自动提交。

整个计算项目存放在根目录\root下，所有MATLAB代码存放在\root\matlab_code下，第i个有限元计算文件存放在\root\abaqus_inp_i下（包括INP文件、UMAT子程序以及提取数据的Python脚本）。

粒子群算法

粒子群算法的实现十分简单，将粒子群规模设定为30，待标定的参数个数为3，函数fitness30用于计算30个粒子的适应度，完整代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
clc;clear;
profile on;%打开探查器, 用统计用时等信息
num=30;%粒子群数
c1=1.5;c2=1.5;%学习因子
N=100;%最大迭代次数

omega0=0.95;omega1=0.4;
omega=@(x) (omega0-omega1)*(N-x)/N+omega1;%线性递减权值

p_limit=[10,30,5;200,100,28];%粒子位置范围
v_limit=[-20,-10,-10;20,10,10];%粒子速度范围

for i=1:num
    p_mat(i,:)=p_limit(1,:)+(p_limit(2,:)-p_limit(1,:)).*rand(1,3);%初始化随机粒子位置
    v_mat(i,:)=v_limit(1,:)+(v_limit(2,:)-v_limit(1,:)).*rand(1,3);%初始化随机粒子速度
end
pbest=p_mat;%粒子的历史最优位置, 即为每个粒子的历史最优位置

%计算初始适应度
fy=fitness30(p_mat);
[min_fy,index]=min(fy);
gbest=pbest(index,:);

fitness_table=zeros(N,1);
%存储初始数据
fid2=fopen('process_data2.txt','at');%初始最优粒子的数据
fprintf(fid2,'%d\t%d\t%.7e\t%.7e\t%.7e\t%.7e\n',...
0,index,min_fy,gbest(1),gbest(2),gbest(3));
fclose(fid2);

%进入迭代
for i=1:N
    
    %更新速度
    for j=1:num
        v_mat(j,:)=omega(i)*v_mat(j,:)+c1*(pbest(j,:)-...
        p_mat(j,:)).*rand(1,3)+c2*(gbest-p_mat(j,:)).*rand(1,3);
        %考虑速度边界
        for k=1:3
            if v_mat(j,k)<v_limit(1,k)
                v_mat(j,k)=v_limit(1,k);
            elseif v_mat(j,k)>v_limit(2,k)
                v_mat(j,k)=v_limit(2,k);
            end
        end
    end
    
    %更新位置
    p_mat=p_mat+v_mat;
    %考虑位置边界
    for j=1:num
        for k=1:3
            if p_mat(j,k)<p_limit(1,k)
                p_mat(j,k)=p_limit(1,k);
            elseif p_mat(j,k)>p_limit(2,k)
                p_mat(j,k)=p_limit(2,k);
            end
        end
    end
    
    fy_new=fitness30(p_mat);%计算新的适应度
    for j=1:num
        if fy_new(j)<fy(j)
            fy(j)=fy_new(j);
            pbest(j,:)=p_mat(j,:);%更新个体最优
        end
    end
    [min_fy,index]=min(fy);
    gbest=pbest(index,:);
    fitness_table(i)=min_fy;%记录每次迭代的全局最优

    %存储数据
    fid2=fopen('process_data2.txt','at');%每次迭代过程最优粒子的数据
    fprintf(fid2,'%d\t%d\t%.7e\t%.7e\t%.7e\t%.7e\n',...
    i,index,min_fy,gbest(1),gbest(2),gbest(3));
    fclose(fid2);

    %第二个终止准则
    if min_fy<0.4
        break;
    end
    
end

profile viewer;

有限元计算

完成有限元计算的三个步骤：前处理（修改参数）、提交计算、后处理（提取数据）

前处理

针对本构参数标定，前处理仅需改变材料参数的数值即可。使用UMAT子程序时，可对ABAQUS的INP文件中进行修改，需要修改的材料参数位于关键字*User Material下。

首先，读取INP文件，找到需要修改位置，并存储在newline中

1
2
3
4
5
6
7
8
9
10
11
12
13
14
absolute_path=[path,'\',InpFile,'.inp'];%生成绝对路径
fid=fopen(absolute_path,'rt+');%以只读模式打开绝对路径下的INP文件
i=0;
newline=cell(1,1000);%预制存储元胞, 长度大于INP文件的行数即可
while ~feof(fid)
    tline=fgetl(fid);
    i=i+1;
    newline{i}=tline;
    if i==408 %遍历到INP文件中需要修改的第408行
        %将该行修改为 1.0,2.0,3.0
        newline{i}=[num2str(1.0),',',num2str(2.0),',',num2str(3.0)];
    end
end
fclose(fid);%关闭绝对路径下的INP文件

接着，根据newline中的数据覆盖原来的INP文件

1
2
3
4
5
fid=fopen(absolute_path,'wt+');%以读写模式打开绝对路径下的INP文件
for k=1:i
    fprintf(fid,'%s\n',newline{k});%覆盖原来的INP文件
end
fclose(fid);%关闭绝对路径下的INP文件

INP文件作为ABAQUS的输入文件，提交后不会被ABAQUS修改，其内容和结构是固定的，因此可在MATLAB中每次遍历修改某些行的内容。

提交计算

ABAQUS提供了提交计算的系统命令，在MATLAB中命令如下

1
2
3
4
5
6
7
8
MatlabPath=pwd();%记下当前Matlab目录, 即\root\matlab_code
cd(path);%进入INP文件所在目录
InpFile = 'Job-1';%INP文件
UserFile = 'my_UMAT';%用户子程序
NumCpu = '1';%计算使用的CPU线程数
InputFile = ['abaqus job=',InpFile,' user=',UserFile,' cpus=',NumCpu];
system(InputFile);%通过系统命令提交ABAQUS计算任务
cd(MatlabPath);%返回Matlab目录, 即\root\matlab_code

后处理

ABAQUS提供了提交计算和后处理脚本的系统命令，在MATLAB中命令如下

1
2
3
4
MatlabPath=pwd();%记下当前Matlab目录, 即\root\matlab_code
cd(path);%进入INP文件所在目录
system('abaqus cae noGUI=get_data.py');%通过系统命令执行Python脚本提取ODB文件中数据
cd(MatlabPath);%返回Matlab目录, 即\root\matlab_code

后处理提出数据还需要编写Python脚本，关键代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
from abaqus import*
from odbAccess import*
import numpy as np
...
odb=openOdb('Job-1.odb')
history=odb.steps['Step-1'].historyRegions
kk=history.keys()
...
for k in kk:
    p=history[k].historyOutputs
...
np.savetxt('data_abaqus.txt',temp_table,fmt='%.6e  %.9e  %.9e')#[t,f,x]
odb.close()

为了提高读取速度，在Python脚本中采用的是历史输出结果historyOutputs，而非场输出结果fieldOutputs，同时需要在ABAQUS中提前设置历史输出量。

计算判断

判断ABAQUS计算是否崩溃、计算提交是否成功、计算是否结束

ABAQUS-Standard求解器允许同时提交多个计算任务，而多任务并行计算时可能会导致某个计算任务崩溃（概率不大）。一旦发生计算崩溃，在粒子群算法中对应的这个粒子也会停止迭代，导致整个计算停止。

在新一轮的\(N\)个计算任务开始前，若文件夹\root\abaqus_inp_i中存在.lck文件，则说明上一轮的计算已经崩溃。

Q：为什么不能是上一轮计算还未完成，导致.lck文件存在？

A：因为一轮的计算正常完成所花费的时间足够完成第i个有限元计算。

为解决该问题，可新建文件夹\root\abaqus_inp_i_OXT（0XT表示新建的第X个的文件夹，前X-1个也都崩溃了），从备份文件夹\root\abaqus_copy中拷贝文件到新建文件夹中，并完成第i个有限元计算，关键代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
%计算失败判断, 更换路径并开启新的计算
if exist([path,'\Job-1.lck'],'file')==2
    temp_str=path;
    if temp_str(end)~='T'%第一次计算失败
        temp_str1=[temp_str,'_01T'];
    else%非第一次计算失败, 最多允许失败九次
        temp_sn=str2double(temp_str(end-2:end-1))+1;
        temp_str1=[temp_str(1:end-3),'0',num2str(temp_sn),'T'];
    end
    mkdir(temp_str1);pause(0.5);%创建新文件夹
    %复制三个文件
    copyfile('\root\abaqus_copy\my_UMAT.for',temp_str1);pause(0.5);
    copyfile('\root\abaqus_copy\get_data.py',temp_str1);pause(0.5);
    copyfile('\root\abaqus_copy\Job-1.inp',temp_str1);pause(0.5);
end

完成计算是否崩溃的判断后，通过系统命令提交ABAQUS进行计算，此时exist([path,'\Job-1.lck'],'file')的值若为2，即提交成功。

当exist([path,'\Job-1.lck'],'file')的值为0时，说明计算结束，可以进行数据提取。

算例结果

计算使用的处理器为\(\mathrm{i}9\)-\(9980\mathrm{XE}\)，具有18核心、36线程，可满足种群数为30的计算需求

种群数为30，最大迭代次数为100，共计完成3000次有限元的前处理、计算以及后处理过程，100轮迭代完成时\(err\approx0.47\)，探查器统计的部分信息如下表所示

函数名称 调用次数 总时间 自用时间
main 1 27008 s 2 s
fitness30 101 27006 s 1724 s
dos 6057 15971 s 15971 s
  • 平均每次过程用时\(\bar{t}=27008~\mathrm{s}/3000 \approx 9~\mathrm{s}\)
  • 该算例迭代到20轮左右时，\(err<1\)
  • 显然，使用粒子群算法来标定材料参数，其标定效果和计算效率都是很可观的