small edits to improve figure 3

Author: GalPerelman

graphs.py

@@ -3,31 +3,52 @@
 
 
 def customize_plot_fossolo(experiment):
-    experiment.plot_constraints = True
+    experiment.plot_constraints = False
     experiment.plot_y_legend_label = ["DeePC"]
     experiment.plot_u_legend_label = ["DeePC"]
     for s in experiment.compare_signals:
         s["plot"] = True
     experiment.plot_legend_cols = 5
     fig = experiment.plot(experiment.wds)
     fig, u, y = experiment.plot_comparison_signals(fig, signals=experiment.compare_signals)
+
+    u, y = experiment.run_pid(kp=0.01, kd=0.1, ki=1, T=experiment.experiment_horizon+experiment.n_train)
+
+    _u = u[-experiment.experiment_horizon:, :]
+    _y = y[-experiment.experiment_horizon:, :]
+
+    mae, me = experiment.get_error(y=_y)
+    ref = np.ones(_y.shape) * experiment.y_ref
+    cost = np.linalg.norm(_y - ref, 'fro') ** 2 + experiment.input_loss * np.linalg.norm(_u, 'fro') ** 2
+    v_count, v_rate = experiment.get_violations(y=_y)
+    print(f"PID --> Cost: {cost:.3f} | MAE: {mae:.3f} | ME: {me:.3f}")
+    print(f"PID --> {experiment.y_lb}-{experiment.y_ub} Violations: {v_count:.0f}"
+          f" | Violations Rate: {v_rate:.3f}")
+    v_count, v_rate = experiment.get_violations(y=_y, lb=20, ub=40)
+    print(f"PID --> {20}-{40} Violations: {v_count:.0f} | Violations Rate: {v_rate:.3f}")
+
     axes = fig.axes
-    axes[0].set_ylim(6, 44)
-    axes[1].set_ylim(26, 65)
+    axes[0].fill_between(range(experiment.n_train, experiment.experiment_horizon+experiment.n_train+2),
+                         experiment.y_lb, experiment.y_ub,
+                         alpha=0.2, color='grey', label="Constraints")
+    axes[0].plot(y, '#35AC78', label="PID", zorder=3, linewidth=1.2)
+    axes[1].step(u, '#35AC78', label="PID", zorder=3, where='post', linewidth=1.2)
+
+    handles, labels = axes[0].get_legend_handles_labels()
+    order = [1, 4, 0, 2, 3, 5]  # reorder the legend items - after adding PID
+    axes[0].legend([handles[i] for i in order], [labels[i] for i in order], ncol=6, columnspacing=0.75, handletextpad=0.2)
+
+    # handles, labels = axes[1].get_legend_handles_labels()
+    # order = [0, 2, 3, 5]  # reorder the legend items - after adding PID
+    # axes[1].legend([handles[i] for i in order], [labels[i] for i in order], ncol=4, columnspacing=0.75, handletextpad=0.15)
+    axes[1].legend(ncol=4, columnspacing=0.75, handletextpad=0.15)
+
+    axes[0].set_ylim(8, 42)
+    axes[1].set_ylim(30, 66)
     axes[2].set_ylim(0.5, 1.6)
     axes[-1].set_xlim(experiment.n_train, experiment.n_train + experiment.experiment_horizon + 4)
     fig.set_size_inches(8, 6)
-    plt.subplots_adjust(left=0.125, right=0.9, bottom=0.1, top=0.9)
-
-    mae, me = experiment.get_error(y=y)
-    ref = np.ones(y.shape) * experiment.y_ref
-    cost = np.linalg.norm(y - ref, 'fro') ** 2 + experiment.input_loss * np.linalg.norm(u, 'fro') ** 2
-    v_count, v_rate = experiment.get_violations(y=y)
-    print(f"{s['name']} --> Cost: {cost:.3f} | MAE: {mae:.3f} | ME: {me:.3f}")
-    print(f"{s['name']} --> {experiment.y_lb}-{experiment.y_ub} Violations: {v_count:.0f}"
-          f" | Violations Rate: {v_rate:.3f}")
-    v_count, v_rate = experiment.get_violations(y=y, lb=20, ub=40)
-    print(f"{s['name']} --> {20}-{40} Violations: {v_count:.0f} | Violations Rate: {v_rate:.3f}")
+    plt.subplots_adjust(left=0.125, right=0.9, bottom=0.1, top=0.9, hspace=0.1)
 
 
 def customize_plot_pescara(experiment):

main.py

@@ -12,8 +12,9 @@
 import graphs
 import system
 import utils
+from pid import PID
 
-COLORS = ["#0077B8", "#DF5353", "#fdc85e", "#b7b3aa"]
+COLORS = ["#0686cc", "#DF5353", "#fdc85e", "#b7b3aa"]#0077B8
 plt.rcParams['axes.prop_cycle'] = plt.cycler(color=COLORS)
 
 
@@ -105,7 +106,7 @@ def run_experiment(self):
 
     def get_error(self, y, n=None):
         if not isinstance(self.y_ref, np.ndarray):
-            y_ref = np.array([self.y_ref])
+            y_ref = np.tile(self.y_ref, (len(y), 1))
         else:
             y_ref = self.y_ref
 
@@ -139,16 +140,16 @@ def plot(self, sys, fig=None, moving_avg_size=0):
             k = 0
 
         if fig is None:
-            inputs_heights = [1 for _ in range(n_inputs + k)]
+            inputs_heights = [1 for _ in range(n_inputs)] + [0.5 for _ in range(k)]
             fig, axes = plt.subplots(nrows=n_inputs + 1 + k, sharex=True, height_ratios=[1] + inputs_heights)
         else:
             axes = fig.axes
 
         t = len(sys.target_values[1:, 0])
         for i, element in enumerate(sys.target_nodes):
-            axes[0].plot(sys.target_values[1:, i], label=self.plot_y_legend_label[i], zorder=4)
+            axes[0].plot(sys.target_values[1:, i], label=self.plot_y_legend_label[i], zorder=10)
 
-        axes[0].hlines(y=self.y_ref, xmin=0, xmax=t, color='k', zorder=5, label="$y_{ref}$")
+        axes[0].hlines(y=self.y_ref, xmin=0, xmax=t, color='k', zorder=1, linestyle='--', label="$y_{ref}$")
         axes[0].axvspan(0, self.n_train, facecolor='grey', alpha=0.3, zorder=0)
         axes[0].grid(True)
         axes[0].set_ylabel(utils.split_label(self.output_y_label, self.plot_y_labels_max_len))
@@ -223,9 +224,20 @@ def plot_comparison_signals(self, fig, signals: list):
             y = sys.target_values[-self.experiment_horizon:, :]
 
             t = len(sys.target_values[1:, 0])
+
+            mae, me = self.get_error(y=y)
+            ref = np.ones(y.shape) * self.y_ref
+            cost = np.linalg.norm(y - ref, 'fro') ** 2 + self.input_loss * np.linalg.norm(u, 'fro') ** 2
+            v_count, v_rate = self.get_violations(y=y)
+            print(f"{s['name']} --> Cost: {cost:.3f} | MAE: {mae:.3f} | ME: {me:.3f}")
+            print(f"{s['name']} --> {self.y_lb}-{self.y_ub} Violations: {v_count:.0f}"
+                  f" | Violations Rate: {v_rate:.3f}")
+            v_count, v_rate = self.get_violations(y=y, lb=20, ub=40)
+            print(f"{s['name']} --> {20}-{40} Violations: {v_count:.0f} | Violations Rate: {v_rate:.3f}")
+
             if "plot" in s and s["plot"]:
                 for i, element in enumerate(sys.target_nodes):
-                    axes[0].plot(sys.target_values[1:, i], label=s["name"], zorder=2)
+                    axes[0].plot(sys.target_values[1:, i], label=s["name"], zorder=2, linewidth=1.2)
                     if self.plot_legend_cols > 0:
                         # axes[0].legend(ncols=self.plot_legend_cols, fontsize=10)
                         handles, labels = axes[0].get_legend_handles_labels()
@@ -234,7 +246,7 @@ def plot_comparison_signals(self, fig, signals: list):
                                        fontsize=9)
 
                 for i, element in enumerate(sys.control_nodes + sys.control_links):
-                    axes[1].step(range(t), sys.implemented[1:, i], label=s["name"], zorder=2, where='post')
+                    axes[1].step(range(t), sys.implemented[1:, i], label=s["name"], zorder=2, where='post', linewidth=1.2)
                     if self.plot_legend_cols > 0:
                         # axes[1].legend(ncols=self.plot_legend_cols, fontsize=10)
                         handles, labels = axes[1].get_legend_handles_labels()
@@ -244,6 +256,29 @@ def plot_comparison_signals(self, fig, signals: list):
 
         return fig, u, y
 
+    def run_pid(self, kp, ki, kd, T):
+        sys = system.WDSControl(inp_path=self.inp_path,
+                                control_links=self.control_links,
+                                control_nodes=self.control_nodes,
+                                target_nodes=self.target_nodes,
+                                target_param=self.target_param
+                                )
+
+        pid = PID(kp=kp, ki=ki, kd=kd, set_point=self.y_ref)
+        current_value = 0.5 * (self.u_lb + self.u_ub)  # Initial value
+        for _ in range(T):
+            u_optimal = pid.compute(current_value=current_value)
+            sys.apply_input(u=np.array([u_optimal]), noise_std=self.noise_std)
+            current_value = sys.get_last_n_samples(1).y[0]
+
+        y = sys.target_values
+        u = sys.implemented
+        ref = np.ones(y.shape) * self.y_ref
+        self.cost = np.linalg.norm(y - ref, 'fro') ** 2 + self.input_loss * np.linalg.norm(u, 'fro') ** 2
+        self.mae, self.me = self.get_error(y=y)
+        self.v_count, self.v_rate = self.get_violations(y=y)
+        return u, y
+
 
 def run_comparable_signal(sys, ref_input_signal, noise_std):
     ref_sys = system.WDSControl(inp_path=sys.inp_path,