How to adjust the parameters to achieve which MV has more action?

I use Q1 and Q2 to control T1, so as to realize the simulation scenario of multi-controls. I want to adjust the parameters to achieve which MV has more action, as shown in the figure. I found that I couldn’t achieve the effect I wanted by adjusting the cost of the MV, can anyone give me some suggestions? thanks!

Here is my code:

import tclab
from tclab import labtime
from tclab import TCLabModel
import numpy as np
import time
import matplotlib.pyplot as plt
from gekko import GEKKO
import json

class tclab_heaterpipe():
    delay_q1_step = 10
    delay_q2_step = 10
    q1_buffer = [0] * delay_q1_step
    q2_buffer = [0] * delay_q2_step
    m = TCLabModel()

    def __init__(self, d1, d2, model):
        if d1 >= 1 and d2 >= 1:
            self.delay_q1_step = int(d1)
            self.delay_q2_step = int(d2)
            self.q1_buffer = [0] * self.delay_q1_step
            self.q2_buffer = [0] * self.delay_q2_step
            self.m = model
        else:
            self.delay_q1_step = 0
            self.delay_q2_step = 0

    def Q1_delay(self, q1):
        if self.delay_q1_step == 0:
            self.m.Q1(q1)
        self.q1_buffer.insert(0, q1)
        self.m.Q1(self.q1_buffer.pop())

    def Q2_delay(self, q2):
        if self.delay_q2_step == 0:
            self.m.Q1(q2)
        self.q2_buffer.insert(0, q2)
        self.m.Q2(self.q2_buffer.pop())

# Connect to Arduino
connected = False
theta = 1
theta2 = 1
T = tclab.setup(connected)
a = T()
tclab_delay = tclab_heaterpipe(theta, theta2, a)
# Turn LED on
print('LED On')
a.LED(100)

# Run time in minutes
run_time = 80.0
# Number of cycles
loops = int(60.0 * run_time)

#########################################################
# Initialize Model
#########################################################
# use remote=True for MacOS
m = GEKKO(name='tclab-mpc', remote=False)

m.time = np.linspace(0, 400, 41)
step = 10
# Temperature (K)
t1sp = 45.0
T1 = np.ones(int(loops / step) + 1) * a.T1  # temperature (degC)
Tsp1 = np.ones(int(loops / step) + 1) * t1sp  # set point (degC)
# heater values
Q1s = np.ones(int(loops / step) + 1) * 0.0
Q2s = np.ones(int(loops / step) + 1) * 0.0

# Parameters
Q1_ss = m.Param(value=0)
TC1_ss = m.Param(value=a.T1)
Q2_ss = m.Param(value=0)
Kp1 = m.Param(value=0.8)
tau1 = m.Param(value=160.0)
Kp2 = m.Param(value=0.1)
tau2 = m.Param(value=160.0)

# Manipulated variable
Q1 = m.MV(value=0, name='q1')
Q1.STATUS = 1  # use to control temperature
Q1.FSTATUS = 0  # no feedback measurement
Q1.LOWER = 0.0
Q1.UPPER = 100.0
Q1.DMAX = 50.0
Q1.DCOST = 1.0
# Q1.COST = 0.1

Q2 = m.MV(value=0, name='q2')
Q2.STATUS = 1  # use to control temperature
Q2.FSTATUS = 0  # no feedback measurement
Q2.LOWER = 0.0
Q2.UPPER = 100.0
Q2.DCOST = 1.0
# Q2.COST = 1.0

# Controlled variable
TC1 = m.CV(value=a.T1, name='tc1')
TC1.STATUS = 1  # minimize error with setpoint range
TC1.FSTATUS = 1  # receive measurement
TC1.TR_INIT = 2  # reference trajectory
TC1.WSPHI = 20
TC1.WSPLO = 20
TC1.TAU = 40  # time constant for response
# TC1.COST = 1

Q1d = m.Var()
m.delay(Q1, Q1d, theta)
Q2d = m.Var()
m.delay(Q2, Q2d, theta2)
# Equation
m.Equation(tau1 * TC1.dt() + (TC1 - TC1_ss) == Kp1 * (Q1d - Q1_ss))
m.Equation(tau2 * TC1.dt() + (TC1 - TC1_ss) == Kp2 * (Q2d - Q2_ss))

# Global Options
m.options.IMODE = 6  # MPC
m.options.CV_TYPE = 3  # Objective type
m.options.NODES = 2  # Collocation nodes
m.options.MAX_TIME = 10
m.options.SOLVER = 1  # 1=APOPT, 3=IPOPT
##################################################################

# Create plot
plt.figure()
plt.ion()
plt.show()

# Main Loop
a.Q1(0)
a.Q2(0)
Q2s[0:] = 0
start_time = time.time()

tm = np.zeros(int(loops / step) + 1)
j = 0

try:
    time_start = time.time()
    labtime_start = labtime.time()
    if (not connected):
        labtime.set_rate(10)
    for i in tclab.clock(loops, adaptive=False):
        i = int(i)
        if (i == 0):
            continue
        print("-----------------------")
        t_real = time.time() - time_start
        t_lab = labtime.time() - labtime_start
        print("real time = {0:4.1f}    lab time = {1:4.1f}    m.time = {1:4.1f}".format(t_real, t_lab, m.time))
        if (i % step != 0):
            continue
        j = i / step
        j = int(j)
        print(j)
        T1[j] = a.T1
        tm[j] = i
        ###############################
        ### MPC CONTROLLER          ###
        ###############################
        TC1.MEAS = T1[j]
        print("T1 meas:{0:4.1f} ".format(a.T1))
        # input setpoint with deadband +/- DT
        DT = 0.5
        TC1.SPHI = Tsp1[j] + DT
        TC1.SPLO = Tsp1[j] - DT

        try:
            # stop model time to solve MPC in cast the solver takes too much time
            if (not connected):
                labtime.stop()
            m.solve(disp=False)
            # start model time
            if (not connected):
                labtime.start()
        except Exception as e:
            if (not connected):
                if (not labtime.running):
                    labtime.start()
            print("sovle's exception:")
            print(e)
            if (j != 0):
                Q1s[j] = Q1s[j - 1]
                Q2s[j] = Q2s[j - 1]
            continue
        # test for successful solution
        if (m.options.APPSTATUS == 1):
            # retrieve the first Q value
            tclab_delay.Q1_delay(Q1.NEWVAL)
            tclab_delay.Q2_delay(Q2.NEWVAL)
            Q1s[j:] = np.ones(len(Q1s) - j) * Q1.NEWVAL
            Q2s[j:] = np.ones(len(Q2s) - j) * Q2.NEWVAL
            # a.Q1(Q1.NEWVAL)
            # a.Q2(Q2.NEWVAL)
            print("Q1 applied with delay: {0:4.1f}  ".format(Q1.NEWVAL))
            print("Q2 applied with delay: {0:4.1f}  ".format(Q2.NEWVAL))
            with open(m.path + '//results.json') as f:
                results = json.load(f)
        else:
            # not successful, set heater to zero
            Q1s[j] = Q1s[j - 1]
            Q2s[j] = Q2s[j - 1]
            print("APPSTATUS is not 1,set Q to 0")

        if (not connected):
            labtime.stop()
        # Plot
        try:
            plt.clf()
            ax = plt.subplot(2, 1, 1)
            ax.grid()
            plt.plot(tm[0:j], T1[0:j], 'ro', markersize=3, label=r'$T_1$')
            plt.plot(tm[0:j], Tsp1[0:j], 'r-', markersize=3, label=r'$T_1 Setpoint$')

            plt.plot(tm[j] + m.time, results['tc1.bcv'], 'r-.', markersize=1, \
                     label=r'$T_1$ predicted', linewidth=1)

            plt.ylabel('Temperature (degC)')
            plt.legend(loc='best')
            ax = plt.subplot(2, 1, 2)
            ax.grid()
            plt.plot(tm[0:j], Q1s[0:j], 'r-', linewidth=3, label=r'$Q_1$')
            plt.plot(tm[0:j], Q2s[0:j], 'b-', linewidth=3, label=r'$Q_2$')
            plt.plot(tm[j] + m.time, Q1.value, 'r-.', \
                     label=r'$Q_1$ plan', linewidth=1)
            plt.plot(tm[j] + m.time, Q2.value, 'b-.', \
                     label=r'$Q_2$ plan', linewidth=1)
            plt.ylabel('Heaters')
            plt.xlabel('Time (sec)')
            plt.legend(loc='best')
            plt.draw()
            plt.pause(0.05)
        except Exception as e:
            print(e)
            pass

        if (not connected):
            labtime.start()

    # Turn off heaters
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    input("Press Enter to continue...")
    a.close()

# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()

# Make sure serial connection still closes when there's an error
except:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Error: Shutting down')
    a.close()
    raise

Solution

There are a many ways to preferentially select one MV over another with additional resources on the course website with Control Objectives and Control Tuning.

COST (minimize or maximize use of MV)
DCOST (change penalty)
DMAX (move change constraint)
TIER (for strong decoupling)

For this particular problem, it appears that the two MVs (Manipulated Variables) Q1 and Q2 values should be a ratio and there is only one CV (Controlled Variable). An easy way to accomplish this is to add another equation such as:

m.Equation(Q2==1.29*Q1)

This consumes a degree a freedom with the equation so that there is effectively only one MV. If the STATUS is off then this could lead to a solver error (Infeasible solution) so you may want to change one of the MVs to a m.Var() type (such as Q2) if this is a problem.