import numpy as np
import matplotlib.pyplot as plt
from APMonitor.apm import *

# write model
model = '''
! apopt MINLP solver options (see apopt.com)
File apopt.opt
 minlp_maximum_iterations 1000     ! minlp iterations
 minlp_max_iter_with_int_sol 50    ! minlp iterations if integer solution is found
 minlp_as_nlp 0                    ! treat minlp as nlp
 nlp_maximum_iterations 200        ! nlp sub-problem max iterations
 minlp_branch_method 1             ! 1 = depth first, 2 = breadth first
 minlp_gap_tol 0.001               ! covergence tolerance
 minlp_integer_tol 0.001           ! maximum deviation from whole number to be considered an integer
 minlp_integer_leaves 0            ! create soft (1) integer leaves or hard (2) integer leaves with branching  
End File

Constants
  c0 = 0.4 
  c1 = 0.2

Parameters
  last

Variables
  x0 = 0.5 , >= 0
  x1 = 0.7 , >= 0
  x2 = 0.0 , >= 0
  int_w = 0 , >= 0 , <= 1

Intermediates
  w = int_w

Equations
  minimize last * x2

  $x0 = x0 - x0*x1 - c0*x0*w 
  $x1 = - x1 + x0*x1 - c1*x1*w
  $x2 = (x0-1)^2 + (x1-1)^2                                                                                       
'''
fid = open('lotka_volterra.apm','w')
fid.write(model)
fid.close()

# write data file
time = np.linspace(0,12,121)
time = np.insert(time, 1, 0.01)
last = np.zeros(122)
last[-1] = 1.0
data = np.vstack((time,last))
np.savetxt('data.csv',data.T,delimiter=',',header='time,last',comments='')

# specify server and application name
s = 'https://byu.apmonitor.com'
#s = 'https://127.0.0.1/'  # for local APMonitor server
a = 'lotka'

apm(s,a,'clear all')
apm_load(s,a,'lotka_volterra.apm')
csv_load(s,a,'data.csv')

apm_option(s,a,'nlc.imode',6)              # Nonlinear control / dynamic optimization
apm_option(s,a,'nlc.nodes',3)

apm_info(s,a,'MV','int_w')                 # M or MV = Manipulated variable - independent variable over time horizon
apm_option(s,a,'int_w.status',1)           # Status: 1=ON, 0=OFF
apm_option(s,a,'int_w.mv_type',0)          # MV Type = Zero Order Hold

apm_option(s,a,'nlc.solver',1)             # 1 = APOPT

# solve
output = apm(s,a,'solve')            
print(output)

# retrieve solution
y = apm_sol(s,a)

plt.figure(1)
plt.step(y['time'],y['int_w'],'r-',label='w (0/1)')
plt.plot(y['time'],y['x0'],'b-',label=r'$x_0$')
plt.plot(y['time'],y['x1'],'k-',label=r'$x_1$')
plt.plot(y['time'],y['x2'],'g-',label=r'$x_2$')
plt.xlabel('Time')
plt.ylabel('Variables')
plt.legend(loc='best')
plt.show()