From ab26dd9be9715456fdeb0504b304b91f72340341 Mon Sep 17 00:00:00 2001 From: pipeacosta Date: Sun, 17 Mar 2024 10:43:42 +0100 Subject: [PATCH] First working version of co-simulation between DPsim and Matlab Signed-off-by: pipeacosta --- .../emt-cosim-udp-matlab-dpsim.py | 141 ++++++++++++++++++ .../emt-cosim-udp-matlab-matlab.m | 94 ++++++++++++ 2 files changed, 235 insertions(+) create mode 100644 examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-dpsim.py create mode 100644 examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-matlab.m diff --git a/examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-dpsim.py b/examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-dpsim.py new file mode 100644 index 0000000000..fbb625b434 --- /dev/null +++ b/examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-dpsim.py @@ -0,0 +1,141 @@ +import sys +import os.path +import logging + +import dpsimpy +import dpsimpyvillas +import numpy as np +import time, socket, errno, signal, struct + +from multiprocessing import Process + +base = os.path.splitext(os.path.basename(sys.argv[0]))[0] +log = logging.getLogger(base) + +if __name__ == '__main__': + local = ("0.0.0.0", 12008) + remote = (os.environ.get('MATLAB_HOST'), 12009) + + logging.basicConfig(format='[%(asctime)s %(name)s %(levelname)s] %(message)s', datefmt='%H:%M:%S', level=logging.DEBUG) + t_0 = 0.0 + t_s = 0.001 + t_f = 1 + + r_1_r = 0.1 + c_1_c = 1 + r_line_r = 0.1 + + sim_name = "EMTCosimUDP1" + + gnd = dpsimpy.emt.SimNode.gnd + n1 = dpsimpy.emt.SimNode("n1") + n2 = dpsimpy.emt.SimNode("n2") + + evs = dpsimpy.emt.ph1.VoltageSource('v_intf', dpsimpy.LogLevel.info) + evs.set_parameters(2) + + r_1 = dpsimpy.emt.ph1.Resistor("r_1", dpsimpy.LogLevel.info) + r_1.set_parameters(r_1_r) + c_1 = dpsimpy.emt.ph1.Capacitor("c_1", dpsimpy.LogLevel.info) + c_1.set_parameters(c_1_c) + r_line = dpsimpy.emt.ph1.Resistor('r_line', dpsimpy.LogLevel.info) + r_line.set_parameters(r_line_r) + + # Initial conditions + n1.set_initial_voltage(5) + n2.set_initial_voltage(2) + + # Connections + r_1.connect([gnd, n1]) + r_line.connect([n2, n1]) + c_1.connect([gnd, n1]) + evs.connect([gnd, n2]) + + sys = dpsimpy.SystemTopology(50, [gnd, n1, n2], [evs, r_1, c_1, r_line]) + + sim = dpsimpy.Simulation(sim_name, loglevel=dpsimpy.LogLevel.debug) + sim.set_domain(dpsimpy.Domain.EMT) + sim.set_system(sys) + sim.set_time_step(t_s) + sim.set_final_time(t_f) + + dpsimpy.Logger.set_log_dir('logs/' + sim_name) + + logger = dpsimpy.Logger(sim_name) + logger.log_attribute('v1', 'v', n1) + logger.log_attribute('v2_S1', 'v', n2) + logger.log_attribute('i', 'i_intf', evs) + logger.log_attribute('ir', 'i_intf', r_line) + logger.log_attribute('V_ref', 'V_ref', evs) + + sim.add_logger(logger) + + n1_v0 = np.array([5.0]) + n2_v0 = np.array([2.0]) + + ir_1_0 = n1_v0 / r_1_r + i_r_line_0 = (n1_v0 - n2_v0) / r_line_r + + r_1.set_intf_voltage(n1_v0) + r_1.set_intf_current(ir_1_0) + c_1.set_intf_voltage(n1_v0) + c_1.set_intf_current(ir_1_0 - i_r_line_0) + r_line.set_intf_voltage(n1_v0 - n2_v0) + r_line.set_intf_current(i_r_line_0) + + evs.set_intf_voltage(n2_v0) + evs.set_intf_current(i_r_line_0) + + sim.start() + + skt = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) + + # This client must be started first + print("Start client...") + + skt.bind(local) + + # Try to connect in case Unix domain socket does not exist yet.. + connected = False + while not connected: + try: + skt.connect(remote) + except socket.error as serr: + if serr.errno not in [errno.ECONNREFUSED, errno.ENOENT]: + raise serr + + print("Not connected. Retrying in 1 sec") + time.sleep(1) + else: + connected = True + + print("Ready. Ctrl-C to quit.") + + + # Gracefully shutdown + def sighandler(signum, frame): + running = False + + + signal.signal(signal.SIGINT, sighandler) + signal.signal(signal.SIGTERM, sighandler) + + running = True + t_k = 0.0 + + while t_k < t_f: + dgram, addr = skt.recvfrom(1024) + if not dgram: + break + else: + u_1 = struct.unpack('f', dgram) + print(u_1) + + sim.get_idobj_attr("v_intf", "V_ref").set(complex(u_1[0],0)) + t_k = sim.next() + y_1 = sim.get_idobj_attr("v_intf", "i_intf").derive_coeff(0,0).get() + skt.send(struct.pack('f', y_1)) + + skt.close() + + print("Bye.") \ No newline at end of file diff --git a/examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-matlab.m b/examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-matlab.m new file mode 100644 index 0000000000..9401950fc3 --- /dev/null +++ b/examples/villas/emt-cosim-udp-matlab/emt-cosim-udp-matlab-matlab.m @@ -0,0 +1,94 @@ +% Author: Andres Acosta +% SPDX-FileCopyrightText: 2014-2024 Institute for Automation of Complex Power Systems, RWTH Aachen University +% SPDX-License-Identifier: Apache-2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +close all; +clc; + +warning('off','Control:analysis:LsimStartTime') +set(0, 'DefaultLineLineWidth', 1.5); + +n_x2 = 1; +n_u2 = 1; +n_y2 = 1; + +R_3 = 1; +C_2 = 1; + +A_2 = -1/(R_3*C_2); +B_2 = 1/C_2; +C_2 = 1; +D_2 = 0; + +x_2_0 = 2; +y_2_0 = C_2*x_2_0; + +fprintf('Output value from S_2: %f\n', y_2_0); + +% Configuration and connection +u = udp(DPSIM_HOST, ... + 'LocalPort', 12009, ... + 'RemotePort', 12008, ... + 'Timeout', 10, ... + 'DatagramTerminateMode', 'on', ... + 'ByteOrder', 'littleEndian' ... +); +disp('UDP Receiver started'); + +fopen(u); +disp('UDP Socket bound'); + +num_values = 1; +t_s = 0.001; +t_f = 1; + +t = 0:t_s:t_f; + +ss_2 = ss(A_2, B_2, C_2, D_2); + +u_2_v = zeros(1, length(t)); +u_2_v(1) = 30; % This should be sent from the other subsystem +y_2_v = zeros(1, length(t)); +y_2_v(1) = y_2_0; + +% Send initial output +fwrite(u, y_2_0, 'float'); + +i = 2; + +while i <= length(t) + % Receive input from S_1 + [ u_2, c_u_2 ] = fread(u, num_values, 'float'); + + fprintf('Received output value from S_1: %f\n', u_2(1)); + + if ~isempty(u_2) + u_2_v(i) = double(u_2(1)); + end + + [y_2,~,x_2] = lsim(ss_2, u_2_v(i-1)*ones(1,2), t(i-1:i), x_2_0); + fprintf('Output value from S_2: %f\n', y_2(end)); + +% fwrite(u, [0, single(y_2(end))], 'float'); + fwrite(u, y_2(end), 'float'); + + x_2_0 = x_2(end); + y_2_v(i) = y_2(end); + i = i + 1; +end + +figure(); +plot(t, u_2_v); +grid on +xlabel('t [s]') +ylabel('$i_2$ [A]', 'Interpreter', 'Latex') + +figure(); +plot(t, y_2_v); +grid on +xlabel('t [s]') +ylabel('$v_2$ [V]', 'Interpreter', 'Latex') + +fclose(u); +delete(u);