Skip to content

Commit 32f2edf

Browse files
bstellatobstellato
committed
Updated also power converter timings
1 parent c2c7bdd commit 32f2edf

File tree

8 files changed

+82
-99
lines changed

8 files changed

+82
-99
lines changed

‎examples/power_converter/power_converter.py

Lines changed: 15 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
from .tail_cost import TailCost
1919
from .quadratic_program import MIQP
2020

21+
2122
class Statistics(object):
2223
def __init__(self, fsw, thd, max_solve_time, min_solve_time,
2324
avg_solve_time, std_solve_time):
@@ -28,6 +29,7 @@ def __init__(self, fsw, thd, max_solve_time, min_solve_time,
2829
self.avg_solve_time = avg_solve_time
2930
self.std_solve_time = std_solve_time
3031

32+
3133
class SimulationResults(object):
3234
"""
3335
Simulation results signals
@@ -42,6 +44,7 @@ def __init__(self, X, U, Y_phase, Y_star_phase, T_e, T_e_des, solve_times):
4244
self.T_e_des = T_e_des
4345
self.solve_times = solve_times
4446

47+
4548
class DynamicalSystem(object):
4649
"""
4750
Power converter dynamical system
@@ -423,7 +426,6 @@ def compute_mpc_input(self, x0, u_prev, solver='gurobi'):
423426

424427
N = qp.N
425428

426-
427429
# Update objective
428430
q = 2. * (qp.q_x.dot(x0) + qp.q_u)
429431

@@ -434,17 +436,17 @@ def compute_mpc_input(self, x0, u_prev, solver='gurobi'):
434436

435437
if solver == 'gurobi':
436438
# Solve problem
437-
prob = mpbpy.QuadprogProblem(qp.P, q, qp.A, qp.l, qp.u, qp.i_idx, x0=u_prev)
439+
prob = mpbpy.QuadprogProblem(qp.P, q, qp.A, qp.l, qp.u, qp.i_idx,
440+
qp.i_l, qp.i_u, x0=u_prev)
438441
res_gurobi = prob.solve(solver=mpbpy.GUROBI, verbose=False,
439442
Threads=1)
440443
u = res_gurobi.x
441444
obj_val = res_gurobi.obj_val
442445
solve_time = res_gurobi.cputime
443446

444-
445447
elif solver == 'miosqp':
446448

447-
if self.solver == None:
449+
if self.solver is None:
448450
# Define problem settings
449451
miosqp_settings = {'eps_int_feas': 1e-02, # integer feasibility tolerance
450452
'max_iter_bb': 2000, # maximum number of iterations
@@ -459,18 +461,14 @@ def compute_mpc_input(self, x0, u_prev, solver='gurobi'):
459461
osqp_settings = {'eps_abs': 1e-03,
460462
'eps_rel': 1e-03,
461463
'eps_prim_inf': 1e-04,
462-
'rho': 0.001,
463-
'sigma': 1e-06,
464-
'alpha': 1.6,
465-
'polish': False,
466-
'max_iter': 2000,
467-
'early_terminate_interval': 25,
464+
# 'rho': 0.001,
465+
'rho': 0.1,
468466
'verbose': False}
469467
self.solver = miosqp.MIOSQP()
470468
self.solver.setup(qp.P, q, qp.A, qp.l,
471-
qp.u, qp.i_idx,
472-
miosqp_settings,
473-
osqp_settings)
469+
qp.u, qp.i_idx, qp.i_l, qp.i_u,
470+
miosqp_settings,
471+
osqp_settings)
474472
else:
475473
self.solver.update_vectors(q, qp.l, qp.u)
476474

@@ -572,7 +570,6 @@ def get_statistics(self, results):
572570
freq_sw = N_sw / (1. / self.params.freq * sim_periods)
573571
fsw = np.mean(freq_sw) # Compute average between 12 switches
574572

575-
576573
# Get THD
577574
t = self.time.t
578575
t_init = init_periods * Nstpp
@@ -590,13 +587,12 @@ def get_statistics(self, results):
590587
max_solve_time, min_solve_time,
591588
avg_solve_time, std_solve_time)
592589

593-
594590
def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
595591
"""
596592
Perform closed loop simulation
597593
"""
598594

599-
print("Simulating closed loop N = %i with solver %s" % \
595+
print("Simulating closed loop N = %i with solver %s" %
600596
(N, solver))
601597

602598
# Reset solver
@@ -613,7 +609,6 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
613609
T_final = self.time.T_final
614610
T_timing = self.time.T_timing
615611

616-
617612
# Compute QP matrices
618613
self.qp_matrices = MIQP(self.dyn_system, N, self.tail_cost)
619614

@@ -633,14 +628,14 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
633628
# Run loop
634629
for i in tqdm(range(T_final)):
635630

636-
637631
# Compute mpc inputs
638632
U[:, i], obj_vals[i], time_temp, u_prev, osqp_iter = \
639-
self.compute_mpc_input(X[:, i], u_prev, solver=solver)
633+
self.compute_mpc_input(X[:, i], u_prev, solver=solver)
640634

641635
# Store time if after the init periods
642636
if i >= self.time.init_periods * self.time.Nstpp:
643-
solve_times[i - self.time.init_periods * self.time.Nstpp] = time_temp
637+
solve_times[i - self.time.init_periods * self.time.Nstpp] = \
638+
time_temp
644639

645640
# Simulate one step
646641
X[:, i + 1], Y[:, i] = self.simulate_one_step(X[:, i], U[:, i])
@@ -656,7 +651,6 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
656651
# Divide total number of average iterations by time steps
657652
miosqp_avg_osqp_iter /= T_final
658653

659-
660654
# Compute additional signals for plotting
661655
Y_phase, Y_star_phase, T_e, T_e_des = self.compute_signals(X)
662656

‎examples/power_converter/quadratic_program.py

Lines changed: 49 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,9 @@
11
import numpy as np
22
import numpy.linalg as nla
3-
import scipy as sp
43
import scipy.linalg as sla
54
import scipy.sparse as spa
65

76

8-
9-
107
class MIQP(object):
118
"""
129
Mixed-Integer Quadratic Program matrices and vectors
@@ -37,7 +34,7 @@ def __init__(self, dyn_system, N, tail_cost):
3734
# A_tilde matrix
3835
A_tilde = np.eye(nx)
3936
for i in range(1, N + 1):
40-
A_tilde = np.vstack((A_tilde, nla.matrix_power(A ,i)))
37+
A_tilde = np.vstack((A_tilde, nla.matrix_power(A, i)))
4138
A_tilde_end = nla.matrix_power(A, N)
4239

4340
# B_tilde matrix
@@ -65,36 +62,32 @@ def __init__(self, dyn_system, N, tail_cost):
6562

6663
Hx = sla.block_diag(Hx, np.zeros((nx, nx)))
6764

68-
6965
# Quadratic form matrices
70-
qp_P = spa.csc_matrix(
71-
2. * (B_tilde.T.dot(Hx).dot(B_tilde) + \
72-
(gamma**N) * B_tilde_end.T.dot(P0).dot(B_tilde_end)))
66+
qp_P = spa.csc_matrix(2. * (B_tilde.T.dot(Hx).dot(B_tilde) +
67+
(gamma**N) *
68+
B_tilde_end.T.dot(P0).dot(B_tilde_end)))
7369
qp_q_x = B_tilde.T.dot(Hx.T).dot(A_tilde) + \
7470
(gamma ** N) * (B_tilde_end).T.dot(P0).dot(A_tilde_end)
7571
qp_q_u = (gamma ** N) * B_tilde_end.T.dot(q0)
7672

7773
# Constant part of the cost function
7874
qp_const_P = A_tilde.T.dot(Hx).dot(A_tilde) + \
79-
(gamma ** N)* A_tilde_end.T.dot(P0).dot(A_tilde_end)
75+
(gamma ** N) * A_tilde_end.T.dot(P0).dot(A_tilde_end)
8076
qp_const_q = ((gamma ** N) * q0.T.dot(A_tilde_end))
8177
qp_const_r = (gamma ** N) * r0
8278

83-
84-
8579
'''
8680
Define constraints
8781
'''
8882

8983
# Matrices required for constraint satisfaction S, R, T
90-
S1 = np.hstack(( np.kron(np.eye(N), W), np.zeros((3 * N, nx))))
91-
S2 = np.hstack(( np.kron(np.eye(N), -W), np.zeros((3 * N, nx))))
92-
S = np.vstack(( S1, S2))
84+
S1 = np.hstack((np.kron(np.eye(N), W), np.zeros((3 * N, nx))))
85+
S2 = np.hstack((np.kron(np.eye(N), -W), np.zeros((3 * N, nx))))
86+
S = np.vstack((S1, S2))
9387

94-
R = np.vstack(( np.kron(np.eye(N), G - T), np.kron(np.eye(N), -G - T)))
88+
R = np.vstack((np.kron(np.eye(N), G - T), np.kron(np.eye(N), -G - T)))
9589
F = np.kron(np.eye(N), T)
9690

97-
9891
# Linear constraints
9992
qp_A = spa.csc_matrix(np.vstack((R - S.dot(B_tilde), F)))
10093

@@ -104,24 +97,23 @@ def __init__(self, dyn_system, N, tail_cost):
10497
# lower bound
10598
qp_l = np.append(-np.inf * np.ones(6 * N), -np.ones(3 * N))
10699

107-
108100
# Constrain bounds to be within -1 and 1
109-
u_sw_idx = np.append(np.ones(3), np.zeros(3))
110-
u_sw_idx = np.tile(u_sw_idx, N)
111-
u_sw_idx = np.flatnonzero(u_sw_idx)
112-
qp_A, qp_l, qp_u = self.add_bounds(u_sw_idx, -1., 1.,
113-
qp_A, qp_l, qp_u)
114-
115-
101+
# u_sw_idx = np.append(np.ones(3), np.zeros(3))
102+
# u_sw_idx = np.tile(u_sw_idx, N)
103+
# u_sw_idx = np.flatnonzero(u_sw_idx)
104+
# qp_A, qp_l, qp_u = self.add_bounds(u_sw_idx, -1., 1.,
105+
# qp_A, qp_l, qp_u)
106+
116107
# SA_tilde needed to update bounds
117108
qp_SA_tilde = S.dot(A_tilde)
118109

119-
120-
121110
# Index of integer variables
122111
i_idx = np.arange(nu * N)
123112

124-
113+
# Bounds on integer variables
114+
# i_idx = u_sw_idx
115+
i_l = -1. * np.ones(nu * N)
116+
i_u = 1. * np.ones(nu * N)
125117

126118
'''
127119
Define problem matrices
@@ -140,32 +132,33 @@ def __init__(self, dyn_system, N, tail_cost):
140132
self.const_r = qp_const_r
141133
self.N = N
142134
self.i_idx = i_idx
143-
144-
145-
def add_bounds(self, i_idx, l_new, u_new, A, l, u):
146-
"""
147-
Add new bounds on the variables
148-
149-
l_new <= x_i <= u_new for i in i_idx
150-
151-
It is done by adding rows to the contraints
152-
153-
l <= A x <= u
154-
"""
155-
156-
n = A.shape[1]
157-
158-
# Enforce integer variables to be binary => {0, 1}
159-
I_int = spa.identity(n).tocsc()
160-
I_int = I_int[i_idx, :]
161-
l_int = np.empty((n,))
162-
l_int.fill(l_new)
163-
l_int = l_int[i_idx]
164-
u_int = np.empty((n,))
165-
u_int.fill(u_new)
166-
u_int = u_int[i_idx]
167-
A = spa.vstack([A, I_int]).tocsc() # Extend problem constraints matrix A
168-
l = np.append(l, l_int) # Extend problem constraints
169-
u = np.append(u, u_int) # Extend problem constraints
170-
171-
return A, l, u
135+
self.i_l = i_l
136+
self.i_u = i_u
137+
138+
# def add_bounds(self, i_idx, l_new, u_new, A, l, u):
139+
# """
140+
# Add new bounds on the variables
141+
#
142+
# l_new <= x_i <= u_new for i in i_idx
143+
#
144+
# It is done by adding rows to the contraints
145+
#
146+
# l <= A x <= u
147+
# """
148+
#
149+
# n = A.shape[1]
150+
#
151+
# # Enforce integer variables to be binary => {0, 1}
152+
# I_int = spa.identity(n).tocsc()
153+
# I_int = I_int[i_idx, :]
154+
# l_int = np.empty((n,))
155+
# l_int.fill(l_new)
156+
# l_int = l_int[i_idx]
157+
# u_int = np.empty((n,))
158+
# u_int.fill(u_new)
159+
# u_int = u_int[i_idx]
160+
# A = spa.vstack([A, I_int]).tocsc() # Extend problem constraints matrix A
161+
# l = np.append(l, l_int) # Extend problem constraints
162+
# u = np.append(u, u_int) # Extend problem constraints
163+
#
164+
# return A, l, u

‎examples/power_converter/run_example.py

Lines changed: 11 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -8,16 +8,14 @@
88
import numpy as np
99
import pandas as pd
1010

11+
# Power converter model files
12+
from .power_converter import Model
1113

1214
# Import plotting library
1315
import matplotlib.pylab as plt
14-
import matplotlib
15-
colors = { 'b': '#1f77b4',
16-
'g': '#2ca02c',
17-
'o': '#ff7f0e'}
18-
19-
# Power converter model files
20-
from .power_converter import Model
16+
colors = {'b': '#1f77b4',
17+
'g': '#2ca02c',
18+
'o': '#ff7f0e'}
2119

2220

2321
def run_example():
@@ -28,16 +26,16 @@ def run_example():
2826
freq = 50. # Switching frequency
2927
torque = 1. # Desired torque
3028
t0 = 0.0 # Initial time
31-
init_periods = 1 # Number of integer period to settle before simulation
29+
init_periods = 1 # Number of periods to settle before simulation
3230
sim_periods = 2 # Numer of simulated periods
3331
flag_steady_trans = 0 # Flag Steady State (0) or Transients (1)
3432

35-
3633
'''
3734
ADP Parameters
3835
'''
3936
gamma = 0.95 # Forgetting factor
4037
N_adp = np.arange(1, 6) # Horizon length
38+
# N_adp = np.array([2])
4139
delta = 5.5 # Switching frequency penalty
4240
N_tail = 50
4341

@@ -46,7 +44,6 @@ def run_example():
4644
k2 = 0.8e03
4745
fsw_des = 300
4846

49-
5047
'''
5148
Setup model
5249
'''
@@ -64,11 +61,10 @@ def run_example():
6461
'''
6562
Allocate output statistics
6663
'''
67-
THD_adp = []
68-
fsw_adp = []
69-
Te_adp = []
70-
Times_adp = []
71-
64+
# THD_adp = []
65+
# fsw_adp = []
66+
# Te_adp = []
67+
# Times_adp = []
7268

7369
'''
7470
Run simulations

‎results/power_converter_currents.pdf

-6 Bytes
Binary file not shown.

‎results/power_converter_inputs.pdf

-584 Bytes
Binary file not shown.

‎results/power_converter_timings.csv

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
,grb_avg,grb_max,grb_min,grb_std,miosqp_avg,miosqp_avg_osqp_iter,miosqp_max,miosqp_min,miosqp_std
2-
1,0.0006160739064216613,0.0038878917694091797,0.00044918060302734375,0.0002451029682178427,0.00018481940031051636,32.66464120370372,0.006685018539428711,5.793571472167969e-05,0.00028991842415354005
3-
2,0.0008126549422740936,0.0025720596313476562,0.0005810260772705078,0.00023366753656506835,0.00019702121615409852,31.08315408549784,0.0013358592987060547,0.00010085105895996094,0.0001669231514913962
4-
3,0.0010482028126716613,0.012631893157958984,0.0006978511810302734,0.000466547408398602,0.00030683636665344236,30.875952519702516,0.0018310546875,0.00013494491577148438,0.0002730323421010312
5-
4,0.0012811090052127839,0.0063550472259521484,0.0008139610290527344,0.0004558778282949365,0.0004734224081039429,34.84607785921957,0.0022630691528320312,0.0001728534698486328,0.00037704791972131347
6-
5,0.0018027685582637787,0.006493091583251953,0.001001119613647461,0.0008631470297120534,0.0006859405338764191,35.02684646162961,0.003576993942260742,0.00021409988403320312,0.0005976953128529478
2+
1,0.0005851049721240998,0.0021419525146484375,0.00043702125549316406,0.00018005306298585615,0.0001693342626094818,25.276909722222214,0.0015380382537841797,8.511543273925781e-05,0.0001441334976546646
3+
2,0.0007725834846496582,0.002134084701538086,0.0005738735198974609,0.00021027299145230218,0.0002366204559803009,25.56510416666666,0.0020699501037597656,0.00010609626770019531,0.00018589408539295232
4+
3,0.0009600065648555755,0.003070831298828125,0.0006809234619140625,0.00031451305184795636,0.00030296996235847475,25.808830171904596,0.002708911895751953,0.00012803077697753906,0.00022741718086372525
5+
4,0.0011823554337024689,0.003644227981567383,0.0007791519165039062,0.0004145728324279434,0.00036506056785583495,26.58755095566407,0.001940011978149414,0.0001518726348876953,0.00026921130796136853
6+
5,0.0016906334459781648,0.006179094314575195,0.0009570121765136719,0.0008155732286066768,0.000550002008676529,30.99985391347512,0.003062009811401367,0.00018286705017089844,0.00042705724358368117

‎results/power_converter_timings.pdf

-18 Bytes
Binary file not shown.

‎run_examples.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,5 +9,5 @@
99

1010
if __name__ == "__main__":
1111

12-
random_miqp.run_example()
13-
# power_converter.run_example()
12+
# random_miqp.run_example()
13+
power_converter.run_example()

0 commit comments

Comments
 (0)