Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import matplotlib.pyplot as plt
- import numpy as np
- from scipy.optimize import least_squares, minimize
- import random
- lambda_ch = True
- sat_u_ch = False
- sat_upri_ch = False
- x_rand = round(random.uniform(0, 1), 1)
- r1 = np.random.rand()*10
- r2 = np.random.rand()*10
- r3 = np.random.rand()*10
- def calculate(args):
- T = 80
- h = 0.01
- iterations = int(T / h)
- t_arr = np.arange(0, T, h)
- x_arr = np.array(t_arr)
- y_arr = np.array(t_arr)
- e_arr = np.array(t_arr)
- u_arr = np.array(t_arr)
- e_sum = 0
- y = 0.1
- x = x_rand
- y_need = 0.3
- e = y_need - x
- u = 0
- sat_arr = [-1, 2]
- for i in range(iterations):
- x_arr[i] = x
- y_arr[i] = y
- e_arr[i] = e
- u_arr[i] = u
- e = y_need - x
- e_sum += + e * h
- e_diff = e - e_arr[i]
- u = args[0] * e + args[1] * e_sum + args[2] * e_diff / h
- if sat_upri_ch:
- u_prev = 0 if i == 0 else u_arr[i]
- change = 0.1
- if abs(u_prev - u) > 0.2:
- u = u_prev + (change if u_prev < u else -1 * change)
- if sat_u_ch:
- if u < sat_arr[0]:
- u = sat_arr[0]
- elif u > sat_arr[1]:
- u = sat_arr[1]
- else:
- pass
- x = x + h * (y)
- y = y + h * (np.sin(x) - u * np.cos(x))
- return e_sum, t_arr, x_arr, y_arr, e_arr, u_arr
- def checkParams(args, max):
- for val in args:
- if abs(val) > max or val < 0:
- return True
- return False
- def optimize(args):
- p_lambda = 1.0
- e_sum, a, b, c, e_arr, e = calculate(args)
- return np.sum(abs(e_arr))+p_lambda*(args[0]**2+args[1]**2+args[2]**2)
- args = minimize(optimize, np.array([r1, r2, r3], dtype=float), method='BFGS') # BFGS
- print(args)
- e_sum, t_arr, x_arr, y_arr, e_arr, u_arr = calculate(args.x)
- plt.subplot(311)
- plt.grid(True)
- plt.plot(t_arr, x_arr, label="y")
- plt.plot(t_arr, y_arr, label="x")
- plt.legend(loc="upper right")
- plt.ylim()
- plt.ylabel("stan")
- plt.subplot(312)
- plt.grid(True)
- plt.plot(t_arr, e_arr)
- plt.ylabel("blad")
- plt.subplot(313)
- plt.grid(True)
- plt.plot(t_arr, u_arr, label="wejscie")
- plt.ylabel("u")
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement