书接上回，《SIMPACK与Python联合仿真——2. C程序代码编写与编译》，继续实现Python代码，并实现联合仿真。

1. Python正弦控制代码编写

将与SIMPACK的通信封装为一个类，类名为SPCKenv，将与SIMPACK的交互抽象为以下类成员函数：

类的初始化：__init__
重置与SIMPACK的联合仿真：reset方法
与SIMPACK的单步交互：step方法

上述类的形式与Gym/Gymnasium的API比较类似，不考虑环境注册等因素是可以与Gym相关代码兼容的。由于SIMPACK本身在Solver->Realtime->Animation中提供了Realtime仿真时的动力学动画效果，因此无需实现用于可视化的render方法。由于设计之初考虑的是多次重复仿真，因此没有像Gym API一样设计close方法，其关闭通信的功能被拆分到step方法的判断结束分支中实现。

SPCKenv类的实现如下：

import time
import numpy as np
import pandas as pd
import socket
import struct

class SPCKenv():
    def __init__(self):
        # 定义最长时间步长
        self.Maxts = 1000 
        # 定义积分步长
        self.CtrlIt = 0.01
        # 计算总仿真时间
        self.SumT = self.Maxts * self.CtrlIt
        # 与C程序之间TCP通信端口
        self.TCP = 9999
        # C程序与SIMPACK之间通信端口
        self.UDP = 12999
    def reset(self):
        # 记录时间步      
        self.nowsteps = 0
        # 以下注释的代码端，用于将SIMPACK启动与Python程序集成至同一Python程序
        # 启动SPCK
        # command = f" /.../Documents/SimpackFiles/99_Blogs/test01 {self.TCP} {self.UDP} {self.CtrlIt} {self.SumT + 10}"
        # pendulum_process = threading.Thread(target=run_command, args=(command,))
        # pendulum_process.start() 
        # 留足给SPCK启动的时间差
        #time.sleep(2)  
        # TCP通信初始化
        SERVER_ADDRESS = '127.0.0.1'  
        self.client_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
        buffer_size = 512 
        self.client_socket.setsockopt(socket.SOL_SOCKET, socket.SO_SNDBUF, buffer_size)
        self.client_socket.setsockopt(socket.SOL_SOCKET, socket.SO_RCVBUF, buffer_size)
        self.client_socket.connect((SERVER_ADDRESS, self.TCP))
        time.sleep(0.5) 
        # SPCK完成与python的TCP连接时需要时间
        OBs = [0, 0, 0, 0]  
        OBs = np.array(OBs)
        # 返回初始观测值
        return OBs, {}
    def step(self,action):
        # 只考虑与SIMPACK联合仿真时，不用关心变量terminated、reward
        terminated = 0 
        truncated = 0
        reward = 0
        self.nowsteps = self.nowsteps + 1
        # 从SIMPACK C程序获取数据
        data = self.client_socket.recv(4)
        ny = struct.unpack('!I', data)[0]
        y_values = []
        for _ in range(ny):
            data = self.client_socket.recv(4)
            y_value = struct.unpack('!f', data)[0]
            y_values.append(y_value)
        # 从SIMPACK中获得6个观测值(y-Outpyuts)
        OBs = [y_values[0],y_values[1], y_values[2], y_values[3], y_values[4], y_values[5]] 
        OBs = np.array(OBs)
        print("Python: 获取SIMPACK的观测值为 ", OBs)   
        if(self.nowsteps >= self.Maxts):
            truncated = 1
        if  terminated or truncated: 
            self.client_socket.close()
            print("Python:关闭socket")
            time.sleep(1.5)     
        else:
            action_value = action.item() if isinstance(action, np.ndarray) else action
            u_values = [action_value, 0, 0, 0]
            u_values_str = " ".join([f"{u_value:.4f}" for u_value in u_values])
            self.client_socket.sendall(f"{u_values_str}\n".encode('utf-8'))
        return OBs, reward, terminated, truncated, {}

代码中，变量terminated、reward与本博客内容无关，无需关心；OBs为从SIMPACK中获得6个观测值(y-Outpyuts)；Python程序通过client_socket.recv从C程序获取数据y-Outpyuts，使用client_socket.sendall将控制力矩u-Inputs发送给SIMPACK。

由于需要对于倒立摆输入时间历程为正弦信号的力矩，因此定义一个根据输入时间t，输出正弦信号的值的Python函数sinInput_signal。正弦的峰值为1，角频率为10，相位差为0。代码如下：

import math
def sinInput_signal(t):
    amplitude = 1  # 峰值
    angular_frequency = 10  # 角频率
    phase_difference = 0  # 相位差
    return amplitude * math.sin(angular_frequency * t + phase_difference)

主程序代码为

Pendulum_Obj = SPCKenv()

# Start a new episode.
obs, info = Pendulum_Obj.reset()
time.sleep(0.5)

# 初始时间
NowTime = 0
# 初始化AllObs列表，用于记录y-Outputs在每次循环中的值
AllObs = []

while True:
    # 输入力矩随时间变化
    action = sinInput_signal(NowTime)
    obs, reward, terminated, truncated, info  = Pendulum_Obj.step(action)
    NowTime = obs[5]
    # 将当前的obs添加到AllObs列表中
    AllObs.append(obs)
    
    if terminated or truncated:
        print("仿真结束\n")
        break;

AllObs变量用于记录在while循环中每次循环的所有观测量obs(y-Outpyuts)，其通过列表来存储所有的obs。每次迭代时，将新的obs添加到列表中。

此外，为了将AllObs中角度、角速度、输入控制力矩随时间历程的变化作图，还定义了下列使用matplotlib的作图代码段落。

from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt

# 将列表转换为ndarray
AllObs = np.array(AllObs)

# 调用Linux的中文字体
font = FontProperties(fname="/usr/share/fonts/truetype/arphic/uming.ttc", size=14)  # 指定字体路径和大小

# 绘制AllObs中第2个变量(旋转角速度)随第5个变量(时间)的变化
plt.figure()
plt.plot(AllObs[:, 5], AllObs[:, 2])
plt.xlabel('Time')
plt.ylabel('theta_dt')
plt.title('旋转角速度随时间的变化', fontproperties=font)
plt.show()

# 绘制AllObs中第3个变量(旋转角度)随第5个变量(时间)的变化
plt.figure()
plt.plot(AllObs[:, 5], AllObs[:, 3])
plt.xlabel('Time')
plt.ylabel('theta')
plt.title('旋转角度随时间的变化', fontproperties=font)
plt.show()

# 绘制AllObs中第4个变量(输入力矩)随第5个变量(时间)的变化
plt.figure()
plt.plot(AllObs[:, 5], AllObs[:, 4])
plt.xlabel('Time')
plt.ylabel('Input Torque')
plt.title('输入力矩随时间的变化', fontproperties=font)
plt.show()

2. 联合仿真流程

首先在编译后的可执行文件test01所在文件夹启动终端1，在TCP端口9999、UDP端口12999上以仿真步长0.01s启动SIMPACK Realtime求解器：

扫描二维码关注公众号，回复： 15295333 查看本文章

./test01 9999 12999 0.01

再在jupyter中运行上述Python主函数，请确保涉及的包在当前Python环境中已被安装。

终端1中将会出现下列提示

username@u2004:...$ ./test01 9999 12999 0.01

Run simpack-rt-slv
PID = 2360390
WARNING: ***********************************************************************************
         * The OLicense server is marked as deprecated in Simpack in 2021x.
         * Only the DS License Server (DSLS) is supported in Simpack 2022 and later.
         * Contact your SIMULIA sales partner to migrate to the DSLS.
         ***********************************************************************************
WARNING: Could not lock memory: Cannot allocate memory.
         Ensure ulimit "max locked memory" of current user is large enough for Simpack Realtime.
Model info: 4 6 0 0 0.000000 /home/yaoyao/Documents/SimpackFiles/99_Blogs/test01.spck.
y(1) = y($Y_cos_theta) = 0.000000
y(2) = y($Y_sin_theta) = -1.000000
y(3) = y($Y_theta_dt) = 0.000000
y(4) = y($Y_theta) = -3.141593
y(5) = y($Y_action) = 0.000000
y(6) = y($Y_runnningtime) = 0.000000
u(1)[1] = $UI_Tq01 (0.000000)
u(2)[2] = $UI_Tq02 (0.000000)
u(3)[3] = $UI_Tq03 (0.000000)
u(4)[4] = $UI_Tq04 (0.000000)
Dimensions of u: 4
Dimensions of y: 6
C:开始calculation loop
C:因Actor客户端主动关闭TCP连接，提前停止Simpack仿真

其中，“C:开始calculation loop”、“C:因Actor客户端主动关闭TCP连接，提前停止Simpack仿真”为上篇博客中C程序自定义的输出，其他为SIMPACK的求解器输出，输出的详细程度受C程序中宏SPCK_VERBOSE的取值控制。

Python代码将会打印各时间步的观测值y-Inputs。

3. 控制效果

通过作图函数，将AllObs各个参量随时间的变化量可视化，如下所示：

或者，在仿真计算过程中，在Solver->Realtime->Animation中开启动画。若求解器正在运行，将通过UDP传输至localhost:11100端口，并在GUI界面中展示。如果仿真速度较快完成，可以在Python代码或C程序中添加延时，以便于观察。当Realtime求解结束，则Realtime Animation也会返回显示原本GUI界面的模型。

SolverSettings中设置如下图。

至此，SIMPACK与Python联合仿真的简单模型就已经全部完成了。

本篇博客为抛砖引玉，后续将会实现更多更复杂的联合仿真器，包括引入真实物理设备（例如方向盘、传感器等），以及将SIMPACK Realtime作为仿真器开发基于数据驱动的控制算法。

SIMPACK与Python联合仿真——3. Python代码与联合仿真

1. Python正弦控制代码编写

2. 联合仿真流程

3. 控制效果

猜你喜欢