Events and simulation time step - continued

A solution to the problem of running certain control algorithms at a faster rate than the simulation is to create a time event list. The closest time event will be the next simulation time instant. This will result in a variable time-step solver and am not sure how this will affect the future network analysis. The user simulation time step therefore becomes the maximum allowable time step. The code is as follows (click on "view raw" below the code box to see the code in a new window):

	# If there has been a system change or time_event
	# or if time is greater than the simulation time step.
	if ((t>=t_ode) or (time_event=="yes")):
	# Update the input values in u matrix
	for c1 in range(len(source_list)):
	component_objects[source_list[c1]].generate_val(source_list, system_loops_copy, sys_mat_e, sys_mat_a, sys_mat_b, sys_mat_u, t, t-t_ode_prev)
	# The ODE solver.
	# Will return the x(k+1) value called
	# as next_state_vec from x(k) value
	# called curr_state_vec
	mat_ode(sys_mat_e, sys_mat_a, sys_mat_b, [curr_state_vec, next_state_vec], sys_mat_u, t-t_ode_prev, ode_var, system_loops_copy, node_list)
	# Update the values of objects based on x(k+1)
	# So far only ammeter affected by this.
	for comps in component_objects.keys():
	component_objects[comps].update_val(system_loop_map, stiff_ratio, sys_mat_e, sys_mat_a, sys_mat_b, next_state_vec, sys_mat_u, branch_params, branch_events)
	# Recalculate the branch currents from the loop currents
	for c1 in range(len(branch_params)):
	branch_params[c1][-1][2]=0.0
	for c2 in range(len(system_loop_map)):
	if system_loop_map[c2][c1]=="forward" or system_loop_map[c2][c1]=="stiff_forward":
	branch_params[c1][-1][2]+=next_state_vec.data[c2][0]
	elif system_loop_map[c2][c1]=="reverse" or system_loop_map[c2][c1]=="stiff_reverse":
	branch_params[c1][-1][2]-=next_state_vec.data[c2][0]
	# x(k)=x(k+1) for next iteration.
	for c1 in range(system_size):
	curr_state_vec.data[c1][0]=next_state_vec.data[c1][0]
	# Save the previous time instant of ODE solver
	t_ode_prev=t
	# If the above code is running because of a
	# change in the system and not because of time
	# greater than t_ode, do not increment the t_ode.
	if time_event=="yes" and t<t_ode:
	time_event="no"
	else:
	t_ode=t+dt
	# Store the measured values.
	if (t>=t_store):
	f.write("%s " %str(t),)
	for c1 in range(len(meter_list)):
	f.write("%s " %component_objects[meter_list[c1]].op_value,)
	#for c1 in range(len(system_loops_copy)):
	# f.write("%s " %next_state_vec.data[c1][0],)
	f.write("\n")
	t_store=t_store+dt_store
	# This time vector will contain all the future time events
	# generated by the control functions and the main ODE solver
	time_vector=[t_ode]
	# Call the control codes only if controlled elements exist.
	if controlled_elements:
	# Call the control functions in the main control programs.
	# Use the eval function to call the functions as string arguments
	for c1 in range(len(control_files)):
	eval("%s(control_file_inputs, control_file_outputs, control_file_staticvars, control_file_timeevents, component_objects, c1, t, time_vector)" %control_functions[c1])
	# The soonest event will be the next time instant.
	time_vector.sort()
	t=time_vector[0]

view raw ode_update.py hosted with ❤ by GitHub

In the control function, one of the last statements updates the time event vector with the time events in the individual control code. So, __control.py will look like this (click on "view raw" below the code box to see the code in a new window):

	import math
	def control_func(interface_inputs, interface_outputs, interface_static, interface_time, circuit_components, pos, t_clock, sys_t_events):
	Vo=circuit_components['9Q'].op_value
	x_tri=interface_static[pos]['x_tri']
	t1=interface_time[pos]['t1']
	import math
	carr_freq=5000.0
	if (x_tri>1.0):
	x_tri=1.0
	if (x_tri<0.0):
	x_tri=0.0
	dt_sample=1.0e-6
	try:
	str(switch_control)
	except:
	switch_control=0.0
	else:
	pass
	if t_clock>=t1:
	x_tri+=(0.5*carr_freq)*dt_sample
	if (x_tri>1.0):
	x_tri=0.0
	if (x_tri<0.5):
	switch_control=1.0
	else:
	switch_control=0.0
	t1=t1+dt_sample
	circuit_components['5C'].control_values[0]=switch_control
	interface_outputs[pos]['5C'][1][2]=switch_control
	interface_static[pos]['x_tri']=x_tri
	interface_time[pos]['t1']=t1
	sys_t_events.append(t1)
	return

view raw control__.py hosted with ❤ by GitHub

The last statement appends the time event t1 to the sys_t_events list. In case they were multiple time events, all of them would be appended to the list.